-
Notifications
You must be signed in to change notification settings - Fork 62
/
coordinate_routines.f90
4267 lines (3930 loc) · 202 KB
/
coordinate_routines.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
!> \file
!> \author Chris Bradley
!> \brief This module contains all coordinate transformation and support routines.
!>
!> \section LICENSE
!>
!> Version: MPL 1.1/GPL 2.0/LGPL 2.1
!>
!> The contents of this file are subject to the Mozilla Public License
!> Version 1.1 (the "License"); you may not use this file except in
!> compliance with the License. You may obtain a copy of the License at
!> http://www.mozilla.org/MPL/
!>
!> Software distributed under the License is distributed on an "AS IS"
!> basis, WITHOUT WARRANTY OF ANY KIND, either express or implied. See the
!> License for the specific language governing rights and limitations
!> under the License.
!>
!> The Original Code is OpenCMISS
!>
!> The Initial Developer of the Original Code is University of Auckland,
!> Auckland, New Zealand, the University of Oxford, Oxford, United
!> Kingdom and King's College, London, United Kingdom. Portions created
!> by the University of Auckland, the University of Oxford and King's
!> College, London are Copyright (C) 2007-2010 by the University of
!> Auckland, the University of Oxford and King's College, London.
!> All Rights Reserved.
!>
!> Contributor(s): Kumar Mithraratne
!>
!> Alternatively, the contents of this file may be used under the terms of
!> either the GNU General Public License Version 2 or later (the "GPL"), or
!> the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
!> in which case the provisions of the GPL or the LGPL are applicable instead
!> of those above. If you wish to allow use of your version of this file only
!> under the terms of either the GPL or the LGPL, and not to allow others to
!> use your version of this file under the terms of the MPL, indicate your
!> decision by deleting the provisions above and replace them with the notice
!> and other provisions required by the GPL or the LGPL. If you do not delete
!> the provisions above, a recipient may use your version of this file under
!> the terms of any one of the MPL, the GPL or the LGPL.
!>
!> This module contains all coordinate transformation and support routines.
MODULE COORDINATE_ROUTINES
USE BaseRoutines
USE Constants
USE CoordinateSystemAccessRoutines
USE INPUT_OUTPUT
USE ISO_VARYING_STRING
USE Kinds
USE Maths
USE Strings
USE Types
#include "macros.h"
IMPLICIT NONE
!!TODO: Should all of the get/set/create/destroy routines be accessed via pointers???
PRIVATE
!Module parameters
!> \addtogroup COORINDATE_ROUTINES_CoordinateSystemTypes COORDINATE_ROUTINES::CoordinateSystemTypes
!> \see COORDINATE_ROUTINES
!> Coordinate system type parameters
!>@{
INTEGER(INTG), PARAMETER :: COORDINATE_RECTANGULAR_CARTESIAN_TYPE=1 !<Rectangular Cartesian coordinate system type \see COORDINATE_ROUTINES_CoordinateSystemTypes,COORDINATE_ROUTINES
INTEGER(INTG), PARAMETER :: COORDINATE_CYLINDRICAL_POLAR_TYPE=2 !<Cylindrical polar coordinate system type \see COORDINATE_ROUTINES_CoordinateSystemTypes,COORDINATE_ROUTINES
INTEGER(INTG), PARAMETER :: COORDINATE_SPHERICAL_POLAR_TYPE=3 !<Spherical polar coordinate system type \see COORDINATE_ROUTINES_CoordinateSystemTypes,COORDINATE_ROUTINES
INTEGER(INTG), PARAMETER :: COORDINATE_PROLATE_SPHEROIDAL_TYPE=4 !<Prolate spheroidal coordinate system type \see COORDINATE_ROUTINES_CoordinateSystemTypes,COORDINATE_ROUTINES
INTEGER(INTG), PARAMETER :: COORDINATE_OBLATE_SPHEROIDAL_TYPE=5 !<Oblate spheroidal coordinate system type \see COORDINATE_ROUTINES_CoordinateSystemTypes,COORDINATE_ROUTINES
!>@}
!> \addtogroup COORDINATE_ROUTINES_RadialInterpolations COORDINATE_ROUTINES::RadialInterpolations
!> \see COORDINATE_ROUTINES
!> \brief The type of radial interpolation for polar coordinate systems
!>@{
INTEGER(INTG), PARAMETER :: COORDINATE_NO_RADIAL_INTERPOLATION_TYPE=0 !<No radial interpolation \see COORDINATE_ROUTINES_RadialInterpolations,COORDINATE_ROUTINES
INTEGER(INTG), PARAMETER :: COORDINATE_RADIAL_INTERPOLATION_TYPE=1 !<r radial interpolation \see COORDINATE_ROUTINES_RadialInterpolations,COORDINATE_ROUTINES
INTEGER(INTG), PARAMETER :: COORDINATE_RADIAL_SQUARED_INTERPOLATION_TYPE=2 !<r^2 radial interpolation \see COORDINATE_ROUTINES_RadialInterpolations,COORDINATE_ROUTINES
INTEGER(INTG), PARAMETER :: COORDINATE_RADIAL_CUBED_INTERPOLATION_TYPE=3 !<r^3 radial interpolation \see COORDINATE_ROUTINES_RadialInterpolations,COORDINATE_ROUTINES
!>@}
!> \addtogroup COORDINATE_ROUTINES_JacobianType COORDINATE_ROUTINES::JacobianType
!> \see COORDINATE_ROUTINES
!> \brief The type of Jacobian to return when coordinate metrics are calculated.
!>@{
INTEGER(INTG), PARAMETER :: COORDINATE_JACOBIAN_NO_TYPE=0 !<No Jacobian \see COORDINATE_ROUTINES_JacobianTypes,COORDINATE_ROUTINES
INTEGER(INTG), PARAMETER :: COORDINATE_JACOBIAN_LINE_TYPE=1 !<Line type Jacobian \see COORDINATE_ROUTINES_JacobianTypes,COORDINATE_ROUTINES
INTEGER(INTG), PARAMETER :: COORDINATE_JACOBIAN_AREA_TYPE=2 !<Area type Jacobian \see COORDINATE_ROUTINES_JacobianTypes,COORDINATE_ROUTINES
INTEGER(INTG), PARAMETER :: COORDINATE_JACOBIAN_VOLUME_TYPE=3 !<Volume type Jacobian \see COORDINATE_ROUTINES_JacobianTypes,COORDINATE_ROUTINES
!>@}
!Module types
!Module variables
CHARACTER(LEN=21) :: COORDINATE_SYSTEM_TYPE_STRING(5) = &
& (/ "Rectangular Cartesian",&
& "Cylindrical Polar ", &
& "Spherical Polar ", &
& "Prolate Spheroidal ", &
& "Oblate Spheroidal " /)
!Interfaces
!>COORDINATE_CONVERT_FROM_RC performs a coordinate transformation from a rectangular cartesian coordinate at the point with
!>coordinate Z(:) to the returned point with coordinate X(:) in the coordinate system identified by COORDINATE_SYSTEM.
INTERFACE COORDINATE_CONVERT_FROM_RC
MODULE PROCEDURE COORDINATE_CONVERT_FROM_RC_DP
MODULE PROCEDURE COORDINATE_CONVERT_FROM_RC_SP
END INTERFACE !COORDINATE_CONVERT_FROM_RC
!>COORDINATE_CONVERT_TO_RC performs a coordinate transformation from a coordinate system identified by COORDINATE_SYSTEM
!>at the point X(:) to the returned point Z(:) in rectangular cartesian coordinates.
INTERFACE COORDINATE_CONVERT_TO_RC
MODULE PROCEDURE COORDINATE_CONVERT_TO_RC_DP
MODULE PROCEDURE COORDINATE_CONVERT_TO_RC_SP
END INTERFACE !COORDINATE_CONVERT_TO_RC
!>Calculates the difference (or delta) between two points in a coordinate system. Discontinuities for polar coordinate
!>systems are accounted for
INTERFACE COORDINATE_DELTA_CALCULATE
MODULE PROCEDURE COORDINATE_DELTA_CALCULATE_DP
!MODULE PROCEDURE COORDINATE_DELTA_CALCULATE_SP
END INTERFACE !COORDINATE_DELTA_CALCULATE
!>Calculates DX(:)/DZ(I) at X, where Z(I) are rectangular cartesian and X(:) are curvilinear coordinates defined by COORDINATE_SYSTEM. \todo CHANGE NAME TO SOMETHING MORE MEANINGFULL?
INTERFACE DXZ
MODULE PROCEDURE DXZ_DP
!MODULE PROCEDURE DXZ_SP
END INTERFACE !Dxz
!!TODO:: CHANGE NAME TO SOMETHING MORE MEANINGFULL?
INTERFACE D2ZX
MODULE PROCEDURE D2ZX_DP
!MODULE PROCEDURE D2ZX_SP
END INTERFACE !D2ZX
!!TODO:: CHANGE NAME TO SOMETHING MORE MEANINGFULL?
INTERFACE DZX
MODULE PROCEDURE DZX_DP
!MODULE PROCEDURE DZX_SP
END INTERFACE !DZX
INTERFACE COORDINATE_DERIVATIVE_CONVERT_TO_RC
MODULE PROCEDURE COORDINATE_DERIVATIVE_CONVERT_TO_RC_DP
MODULE PROCEDURE COORDINATE_DERIVATIVE_CONVERT_TO_RC_SP
END INTERFACE !COORDINATE_DERIVATIVE_CONVERT_TO_RC
PUBLIC COORDINATE_RECTANGULAR_CARTESIAN_TYPE,COORDINATE_CYLINDRICAL_POLAR_TYPE,COORDINATE_SPHERICAL_POLAR_TYPE, &
& COORDINATE_PROLATE_SPHEROIDAL_TYPE,COORDINATE_OBLATE_SPHEROIDAL_TYPE
PUBLIC COORDINATE_NO_RADIAL_INTERPOLATION_TYPE,COORDINATE_RADIAL_INTERPOLATION_TYPE, &
& COORDINATE_RADIAL_SQUARED_INTERPOLATION_TYPE,COORDINATE_RADIAL_CUBED_INTERPOLATION_TYPE
PUBLIC COORDINATE_JACOBIAN_NO_TYPE,COORDINATE_JACOBIAN_LINE_TYPE,COORDINATE_JACOBIAN_AREA_TYPE,COORDINATE_JACOBIAN_VOLUME_TYPE
PUBLIC COORDINATE_SYSTEM_TYPE_STRING
PUBLIC COORDINATE_CONVERT_FROM_RC,COORDINATE_CONVERT_TO_RC,COORDINATE_DELTA_CALCULATE,COORDINATE_DERIVATIVE_NORM, &
& COORDINATE_INTERPOLATION_ADJUST,COORDINATE_INTERPOLATION_PARAMETERS_ADJUST,COORDINATE_METRICS_CALCULATE
PUBLIC COORDINATE_SYSTEM_DIMENSION_GET,COORDINATE_SYSTEM_DIMENSION_SET
PUBLIC COORDINATE_SYSTEM_FOCUS_GET,COORDINATE_SYSTEM_FOCUS_SET
PUBLIC Coordinates_RadialInterpolationTypeGet,Coordinates_RadialInterpolationTypeSet
PUBLIC COORDINATE_SYSTEM_TYPE_GET,COORDINATE_SYSTEM_TYPE_SET
PUBLIC COORDINATE_SYSTEM_ORIGIN_GET,COORDINATE_SYSTEM_ORIGIN_SET
PUBLIC COORDINATE_SYSTEM_ORIENTATION_GET,COORDINATE_SYSTEM_ORIENTATION_SET
PUBLIC COORDINATE_SYSTEM_CREATE_START,COORDINATE_SYSTEM_CREATE_FINISH
PUBLIC COORDINATE_SYSTEM_DESTROY
PUBLIC COORDINATE_DERIVATIVE_CONVERT_TO_RC
PUBLIC COORDINATE_SYSTEMS_INITIALISE,COORDINATE_SYSTEMS_FINALISE
PUBLIC Coordinates_MaterialSystemCalculate
CONTAINS
!
!================================================================================================================================
!
!>Performs a coordinate transformation from a rectangular cartesian coordinate at the point with coordinate Z(:) to the returned point with coordinate X(:) in the coordinate system identified by COORDINATE_SYSTEM for double precision coordinates.
FUNCTION COORDINATE_CONVERT_FROM_RC_DP(COORDINATE_SYSTEM,Z,ERR,ERROR)
!Argument variables
TYPE(COORDINATE_SYSTEM_TYPE), INTENT(IN) :: COORDINATE_SYSTEM !<The coordinate system to perform the conversion on
REAL(DP), INTENT(IN) :: Z(:) !<The rectangular cartesian coordiantes to convert
INTEGER(INTG), INTENT(OUT) :: ERR !<The error code
TYPE(VARYING_STRING), INTENT(OUT) :: ERROR !<The error string
!Function variable
REAL(DP) :: COORDINATE_CONVERT_FROM_RC_DP(SIZE(Z,1))
!Local variables
REAL(DP) :: A1,A2,A3,A4,A5,A6,A7,A8,A9,FOCUS
ENTERS("COORDINATE_CONVERT_FROM_RC_DP",ERR,ERROR,*999)
COORDINATE_CONVERT_FROM_RC_DP=0.0_DP
IF(SIZE(Z,1)<COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS) &
& CALL FlagError("Size of Z is less than the number of dimensions.",ERR,ERROR,*999)
SELECT CASE(COORDINATE_SYSTEM%TYPE)
CASE(COORDINATE_RECTANGULAR_CARTESIAN_TYPE)
COORDINATE_CONVERT_FROM_RC_DP(1:COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS)=Z(1:COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS)
CASE(COORDINATE_CYLINDRICAL_POLAR_TYPE)
SELECT CASE(COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS)
CASE(2)
COORDINATE_CONVERT_FROM_RC_DP(1)=SQRT(Z(1)**2+Z(2)**2)
COORDINATE_CONVERT_FROM_RC_DP(2)=ATAN2(Z(1),Z(2))
CASE(3)
COORDINATE_CONVERT_FROM_RC_DP(1)=SQRT(Z(1)**2+Z(2)**2)
COORDINATE_CONVERT_FROM_RC_DP(2)=ATAN2(Z(1),Z(2))
COORDINATE_CONVERT_FROM_RC_DP(3)=Z(3)
CASE DEFAULT
CALL FlagError("Invalid number of coordinates.",ERR,ERROR,*999)
END SELECT
IF(COORDINATE_CONVERT_FROM_RC_DP(2)<0.0_DP) &
& COORDINATE_CONVERT_FROM_RC_DP(2)=COORDINATE_CONVERT_FROM_RC_DP(2)+2.0_DP*PI !reference coordinate 0->2*pi
CASE(COORDINATE_SPHERICAL_POLAR_TYPE)
IF(COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS==3) THEN
COORDINATE_CONVERT_FROM_RC_DP(1)=SQRT(Z(1)**2+Z(2)**2+Z(3)**2)
IF(ABS(Z(1))>=ZERO_TOLERANCE.OR.ABS(Z(2))>=ZERO_TOLERANCE) THEN
COORDINATE_CONVERT_FROM_RC_DP(2)=ATAN2(Z(2),Z(1))
ELSE
COORDINATE_CONVERT_FROM_RC_DP(2)=0.0_DP
ENDIF
A1=SQRT(Z(1)**2+Z(2)**2)
IF(ABS(Z(3))>=ZERO_TOLERANCE.OR.ABS(A1)>=ZERO_TOLERANCE) THEN
COORDINATE_CONVERT_FROM_RC_DP(3)=ATAN2(Z(3),A1)
ELSE
COORDINATE_CONVERT_FROM_RC_DP(3)=0.0_DP
ENDIF
ELSE
CALL FlagError("Invalid number of coordinates.",ERR,ERROR,*999)
ENDIF
CASE(COORDINATE_PROLATE_SPHEROIDAL_TYPE)
FOCUS=COORDINATE_SYSTEM%FOCUS
IF(COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS==3) THEN
A1=Z(1)**2+Z(2)**2+Z(3)**2-FOCUS**2
A2=SQRT(A1**2+4.0_DP*(FOCUS**2)*(Z(2)**2+Z(3)**2))
A3=2.0_DP*FOCUS**2
A4=MAX((A2+A1)/A3,0.0_DP)
A5=MAX((A2-A1)/A3,0.0_DP)
A6=SQRT(A4)
A7=MIN(SQRT(A5),1.0_DP)
IF(ABS(A7)<=1.0_DP) THEN
A8=ASIN(A7)
ELSE
A8=0.0_DP
CALL FLAG_WARNING("Put A8=0 since ABS(A8)>1.",ERR,ERROR,*999)
ENDIF
IF((ABS(Z(3))<ZERO_TOLERANCE).OR.(ABS(A6)<ZERO_TOLERANCE).OR.(ABS(A7)<ZERO_TOLERANCE)) THEN
A9=0.0_DP
ELSE
IF(ABS(A6*A7)>0.0_DP) THEN
A9=Z(3)/(FOCUS*A6*A7)
ELSE
A9=0.0_DP
CALL FLAG_WARNING("Put A9=0 since A6*A7=0.",ERR,ERROR,*999)
ENDIF
IF(A9>=1.0_DP) THEN
A9=PI/2.0_DP
ELSE IF(A9<=-1.0_DP) THEN
A9=-PI/2.0_DP
ELSE
A9=ASIN(A9)
ENDIF
ENDIF
COORDINATE_CONVERT_FROM_RC_DP(1)=LOG(A6+SQRT(A4+1.0_DP))
IF(ABS(Z(1))>=ZERO_TOLERANCE) THEN
COORDINATE_CONVERT_FROM_RC_DP(2)=A8
ELSE
COORDINATE_CONVERT_FROM_RC_DP(2)=PI-A8
ENDIF
IF(ABS(Z(2))>ZERO_TOLERANCE) THEN
COORDINATE_CONVERT_FROM_RC_DP(3)=MOD(A9+2.0_DP*PI,2.0_DP*PI)
ELSE
COORDINATE_CONVERT_FROM_RC_DP(3)=PI-A9
ENDIF
ELSE
CALL FlagError("Invalid number of coordinates.",ERR,ERROR,*999)
ENDIF
CASE(COORDINATE_OBLATE_SPHEROIDAL_TYPE)
CALL FlagError("Not implemented.",ERR,ERROR,*999)
CASE DEFAULT
CALL FlagError("Invalid coordinate type.",ERR,ERROR,*999)
END SELECT
EXITS("COORDINATE_CONVERT_FROM_RC_DP")
RETURN
999 ERRORSEXITS("COORDINATE_CONVERT_FROM_RC_DP",ERR,ERROR)
RETURN
END FUNCTION COORDINATE_CONVERT_FROM_RC_DP
!
!================================================================================================================================
!
!>COORDINATE_CONVERT_FROM_RC_SP performs a coordinate transformation from a rectangular cartesian coordinate at the
!>point with coordinate Z(:) to the returned point with coordinate X(:) in the coordinate system identified by
!>COORDINATE_SYSTEM for single precision coordinates.
FUNCTION COORDINATE_CONVERT_FROM_RC_SP(COORDINATE_SYSTEM,Z,ERR,ERROR)
!Argument variables
TYPE(COORDINATE_SYSTEM_TYPE), INTENT(IN) :: COORDINATE_SYSTEM !<The coordinate system to convert from RC to
REAL(SP), INTENT(IN) :: Z(:) !<The coordinate to convert from rectangular cartesian
INTEGER(INTG), INTENT(OUT) :: ERR !<The error code
TYPE(VARYING_STRING), INTENT(OUT) :: ERROR !<The error string
!Function variable
REAL(SP) :: COORDINATE_CONVERT_FROM_RC_SP(SIZE(Z,1))
!Local variables
REAL(SP) :: A1,A2,A3,A4,A5,A6,A7,A8,A9,FOCUS
ENTERS("COORDINATE_CONVERT_FROM_RC_SP",ERR,ERROR,*999)
COORDINATE_CONVERT_FROM_RC_SP=0.0_SP
IF(SIZE(Z,1)<COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS) &
& CALL FlagError("Size of Z is less than the number of dimensions.",ERR,ERROR,*999)
SELECT CASE(COORDINATE_SYSTEM%TYPE)
CASE(COORDINATE_RECTANGULAR_CARTESIAN_TYPE)
COORDINATE_CONVERT_FROM_RC_SP(1:COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS)=Z(1:COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS)
CASE(COORDINATE_CYLINDRICAL_POLAR_TYPE)
SELECT CASE(COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS)
CASE(2)
COORDINATE_CONVERT_FROM_RC_SP(1)=SQRT(Z(1)**2+Z(2)**2)
COORDINATE_CONVERT_FROM_RC_SP(2)=ATAN2(Z(1),Z(2))
CASE(3)
COORDINATE_CONVERT_FROM_RC_SP(1)=SQRT(Z(1)**2+Z(2)**2)
COORDINATE_CONVERT_FROM_RC_SP(2)=ATAN2(Z(1),Z(2))
COORDINATE_CONVERT_FROM_RC_SP(3)=Z(3)
CASE DEFAULT
CALL FlagError("Invalid number of coordinates.",ERR,ERROR,*999)
END SELECT
IF(COORDINATE_CONVERT_FROM_RC_SP(2)<0.0_SP) &
& COORDINATE_CONVERT_FROM_RC_SP(2)=COORDINATE_CONVERT_FROM_RC_SP(2)+2.0_SP*REAL(PI,SP) !reference coordinate 0->2*pi
CASE(COORDINATE_SPHERICAL_POLAR_TYPE)
IF(COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS==3) THEN
COORDINATE_CONVERT_FROM_RC_SP(1)=SQRT(Z(1)**2+Z(2)**2+Z(3)**2)
IF(ABS(Z(1))>=ZERO_TOLERANCE_SP.OR.ABS(Z(2))>ZERO_TOLERANCE_SP) THEN
COORDINATE_CONVERT_FROM_RC_SP(2)=ATAN2(Z(2),Z(1))
ELSE
COORDINATE_CONVERT_FROM_RC_SP(2)=0.0_SP
ENDIF
A1=SQRT(Z(1)**2+Z(2)**2)
IF(ABS(Z(3))>=ZERO_TOLERANCE_SP.OR.ABS(A1)>ZERO_TOLERANCE_SP) THEN
COORDINATE_CONVERT_FROM_RC_SP(3)=ATAN2(Z(3),A1)
ELSE
COORDINATE_CONVERT_FROM_RC_SP(3)=0.0_SP
ENDIF
ELSE
CALL FlagError("Invalid number of coordinates.",ERR,ERROR,*999)
ENDIF
CASE(COORDINATE_PROLATE_SPHEROIDAL_TYPE)
FOCUS=REAL(COORDINATE_SYSTEM%FOCUS,SP)
IF(COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS==3) THEN
A1=Z(1)**2+Z(2)**2+Z(3)**2-FOCUS**2
A2=SQRT(A1**2+4.0_SP*(FOCUS**2)*(Z(2)**2+Z(3)**2))
A3=2.0_SP*FOCUS**2
A4=MAX((A2+A1)/A3,0.0_SP)
A5=MAX((A2-A1)/A3,0.0_SP)
A6=SQRT(A4)
A7=MIN(SQRT(A5),1.0_SP)
IF(ABS(A7)<=1.0_SP) THEN
A8=ASIN(A7)
ELSE
A8=0.0_SP
CALL FLAG_WARNING("Put A8=0 since ABS(A8)>1.",ERR,ERROR,*999)
ENDIF
IF((ABS(Z(3))<ZERO_TOLERANCE_SP).OR.(ABS(A6)<ZERO_TOLERANCE_SP).OR.(ABS(A7)<ZERO_TOLERANCE_SP)) THEN
A9=0.0_SP
ELSE
IF(ABS(A6*A7)>ZERO_TOLERANCE_SP) THEN
A9=Z(3)/(FOCUS*A6*A7)
ELSE
A9=0.0_SP
CALL FLAG_WARNING("Put A9=0 since A6*A7=0.",ERR,ERROR,*999)
ENDIF
IF(A9>=1.0_SP) THEN
A9=REAL(PI,SP)/2.0_SP
ELSE IF(A9<=-1.0_SP) THEN
A9=-REAL(PI,SP)/2.0_SP
ELSE
A9=ASIN(A9)
ENDIF
ENDIF
COORDINATE_CONVERT_FROM_RC_SP(1)=LOG(A6+SQRT(A4+1.0_SP))
IF(ABS(Z(1))>=ZERO_TOLERANCE_SP) THEN
COORDINATE_CONVERT_FROM_RC_SP(2)=A8
ELSE
COORDINATE_CONVERT_FROM_RC_SP(2)=REAL(PI,SP)-A8
ENDIF
IF(ABS(Z(2))>=ZERO_TOLERANCE_SP) THEN
COORDINATE_CONVERT_FROM_RC_SP(3)=MOD(A9+2.0_SP*REAL(PI,SP),2.0_SP*&
& REAL(PI,SP))
ELSE
COORDINATE_CONVERT_FROM_RC_SP(3)=REAL(PI,SP)-A9
ENDIF
ELSE
CALL FlagError("Invalid number of coordinates.",ERR,ERROR,*999)
ENDIF
CASE(COORDINATE_OBLATE_SPHEROIDAL_TYPE)
CALL FlagError("Not implemented.",ERR,ERROR,*999)
CASE DEFAULT
CALL FlagError("Invalid coordinate type.",ERR,ERROR,*999)
END SELECT
EXITS("COORDINATE_CONVERT_FROM_RC_SP")
RETURN
999 ERRORSEXITS("COORDINATE_CONVERT_FROM_RC_SP",ERR,ERROR)
RETURN
END FUNCTION COORDINATE_CONVERT_FROM_RC_SP
!
!================================================================================================================================
!
!>COORDINATE_CONVERT_TO_RC_DP performs a coordinate transformation from a coordinate system identified by
!>COORDINATE_SYSTEM at the point X(:) to the returned point Z(:) in rectangular cartesian coordinates for
!>double precision coordinates.
FUNCTION COORDINATE_CONVERT_TO_RC_DP(COORDINATE_SYSTEM,X,ERR,ERROR)
!Argument variables
TYPE(COORDINATE_SYSTEM_TYPE), INTENT(IN) :: COORDINATE_SYSTEM !<The coordinate system to convert to rectangular cartesian
REAL(DP), INTENT(IN) :: X(:) !<The coordiante to convert
INTEGER(INTG), INTENT(OUT) :: ERR !<The error coode
TYPE(VARYING_STRING), INTENT(OUT) :: ERROR !<The error string
!Function variable
REAL(DP) :: COORDINATE_CONVERT_TO_RC_DP(SIZE(X,1))
!Local variables
REAL(DP) :: FOCUS
ENTERS("COORDINATE_CONVERT_TO_RC_DP",ERR,ERROR,*999)
COORDINATE_CONVERT_TO_RC_DP=0.0_DP
IF(SIZE(X,1)<COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS) &
& CALL FlagError("Size of X is less than the number of dimensions.",ERR,ERROR,*999)
SELECT CASE(COORDINATE_SYSTEM%TYPE)
CASE(COORDINATE_RECTANGULAR_CARTESIAN_TYPE)
COORDINATE_CONVERT_TO_RC_DP(1:COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS)=X(1:COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS)
CASE(COORDINATE_CYLINDRICAL_POLAR_TYPE)
SELECT CASE(COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS)
CASE(2)
COORDINATE_CONVERT_TO_RC_DP(1)=X(1)*COS(X(2))
COORDINATE_CONVERT_TO_RC_DP(2)=X(1)*SIN(X(2))
CASE(3)
COORDINATE_CONVERT_TO_RC_DP(1)=X(1)*COS(X(2))
COORDINATE_CONVERT_TO_RC_DP(2)=X(1)*SIN(X(2))
COORDINATE_CONVERT_TO_RC_DP(3)=X(3)
CASE DEFAULT
CALL FlagError("Invalid number of coordinates.",ERR,ERROR,*999)
END SELECT
CASE(COORDINATE_SPHERICAL_POLAR_TYPE)
IF(COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS==3) THEN
COORDINATE_CONVERT_TO_RC_DP(1)=X(1)*COS(X(2))*COS(X(3))
COORDINATE_CONVERT_TO_RC_DP(2)=X(1)*SIN(X(2))*COS(X(3))
COORDINATE_CONVERT_TO_RC_DP(3)=X(1)*SIN(X(3))
ELSE
CALL FlagError("Invalid number of coordinates.",ERR,ERROR,*999)
ENDIF
CASE(COORDINATE_PROLATE_SPHEROIDAL_TYPE)
IF(COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS==3) THEN
FOCUS=COORDINATE_SYSTEM%FOCUS
COORDINATE_CONVERT_TO_RC_DP(1)=FOCUS*COSH(X(1))*COS(X(2))
COORDINATE_CONVERT_TO_RC_DP(2)=FOCUS*SINH(X(1))*SIN(X(2))*COS(X(3))
COORDINATE_CONVERT_TO_RC_DP(3)=FOCUS*SINH(X(1))*SIN(X(2))*SIN(X(3))
ELSE
CALL FlagError("Invalid number of coordinates.",ERR,ERROR,*999)
ENDIF
CASE(COORDINATE_OBLATE_SPHEROIDAL_TYPE)
IF(COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS==3) THEN
FOCUS=COORDINATE_SYSTEM%FOCUS
COORDINATE_CONVERT_TO_RC_DP(1)=FOCUS*COSH(X(1))*COS(X(2))*COS(X(3))
COORDINATE_CONVERT_TO_RC_DP(2)=FOCUS*SINH(X(1))*SIN(X(2))
COORDINATE_CONVERT_TO_RC_DP(3)=FOCUS*COSH(X(1))*COS(X(2))*SIN(X(3))
ELSE
CALL FlagError("Invalid number of coordinates.",ERR,ERROR,*999)
ENDIF
CASE DEFAULT
CALL FlagError("Invalid coordinate type.",ERR,ERROR,*999)
END SELECT
EXITS("COORDINATE_CONVERT_TO_RC_DP")
RETURN
999 ERRORSEXITS("COORDINATE_CONVERT_TO_RC_DP",ERR,ERROR)
RETURN
END FUNCTION COORDINATE_CONVERT_TO_RC_DP
!
!================================================================================================================================
!
!>COORDINATE_CONVERT_TO_RC_SP performs a coordinate transformation from a coordinate system identified by
!>COORDINATE_SYSTEM at the point X(:) to the returned point Z(:) in rectangular cartesian coordinates for
!>single precision coordinates.
FUNCTION COORDINATE_CONVERT_TO_RC_SP(COORDINATE_SYSTEM,X,ERR,ERROR)
!Argument variables
TYPE(COORDINATE_SYSTEM_TYPE), INTENT(IN) :: COORDINATE_SYSTEM !<The coordinate system to convert to rectangular cartesian
REAL(SP), INTENT(IN) :: X(:) !<The coordinate to convert
INTEGER(INTG), INTENT(OUT) :: ERR !<The error code
TYPE(VARYING_STRING), INTENT(OUT) :: ERROR !<The error string
!Function variable
REAL(SP) :: COORDINATE_CONVERT_TO_RC_SP(SIZE(X,1))
!Local variables
REAL(SP) :: FOCUS
ENTERS("COORDINATE_CONVERT_TO_RC_SP",ERR,ERROR,*999)
COORDINATE_CONVERT_TO_RC_SP=0.0_SP
IF(SIZE(X,1)<COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS) &
& CALL FlagError("Size of X is less than the number of dimensions.",ERR,ERROR,*999)
SELECT CASE(COORDINATE_SYSTEM%TYPE)
CASE(COORDINATE_RECTANGULAR_CARTESIAN_TYPE)
COORDINATE_CONVERT_TO_RC_SP(1:COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS)=X(1:COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS)
CASE(COORDINATE_CYLINDRICAL_POLAR_TYPE)
SELECT CASE(COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS)
CASE(2)
COORDINATE_CONVERT_TO_RC_SP(1)=X(1)*COS(X(2))
COORDINATE_CONVERT_TO_RC_SP(2)=X(1)*SIN(X(2))
CASE(3)
COORDINATE_CONVERT_TO_RC_SP(1)=X(1)*COS(X(2))
COORDINATE_CONVERT_TO_RC_SP(2)=X(1)*SIN(X(2))
COORDINATE_CONVERT_TO_RC_SP(3)=X(3)
CASE DEFAULT
CALL FlagError("Invalid number of coordinates.",ERR,ERROR,*999)
END SELECT
CASE(COORDINATE_SPHERICAL_POLAR_TYPE)
IF(COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS==3) THEN
COORDINATE_CONVERT_TO_RC_SP(1)=X(1)*COS(X(2))*COS(X(3))
COORDINATE_CONVERT_TO_RC_SP(2)=X(1)*SIN(X(2))*COS(X(3))
COORDINATE_CONVERT_TO_RC_SP(3)=X(1)*SIN(X(3))
ELSE
CALL FlagError("Invalid number of coordinates.",ERR,ERROR,*999)
ENDIF
CASE(COORDINATE_PROLATE_SPHEROIDAL_TYPE)
IF(COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS==3) THEN
FOCUS=REAL(COORDINATE_SYSTEM%FOCUS,SP)
COORDINATE_CONVERT_TO_RC_SP(1)=FOCUS*COSH(X(1))*COS(X(2))
COORDINATE_CONVERT_TO_RC_SP(2)=FOCUS*SINH(X(1))*SIN(X(2))*COS(X(3))
COORDINATE_CONVERT_TO_RC_SP(3)=FOCUS*SINH(X(1))*SIN(X(2))*SIN(X(3))
ELSE
CALL FlagError("Invalid number of coordinates.",ERR,ERROR,*999)
ENDIF
CASE(COORDINATE_OBLATE_SPHEROIDAL_TYPE)
IF(COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS==3) THEN
FOCUS=REAL(COORDINATE_SYSTEM%FOCUS,SP)
COORDINATE_CONVERT_TO_RC_SP(1)=FOCUS*COSH(X(1))*COS(X(2))*COS(X(3))
COORDINATE_CONVERT_TO_RC_SP(2)=FOCUS*SINH(X(1))*SIN(X(2))
COORDINATE_CONVERT_TO_RC_SP(3)=FOCUS*COSH(X(1))*COS(X(2))*SIN(X(3))
ELSE
CALL FlagError("Invalid number of coordinates.",ERR,ERROR,*999)
ENDIF
CASE DEFAULT
CALL FlagError("Invalid coordinate type.",ERR,ERROR,*999)
END SELECT
EXITS("COORDINATE_CONVERT_TO_RC_SP")
RETURN
999 ERRORSEXITS("COORDINATE_CONVERT_TO_RC_SP",ERR,ERROR)
RETURN
END FUNCTION COORDINATE_CONVERT_TO_RC_SP
!
!================================================================================================================================
!
!>Calculates the difference (or detlta) between the point X and the point Y i.e., Y-X, in the given coordinate system.
!>0->2Pi discontinuities with polar coordinates are accounted for.
FUNCTION COORDINATE_DELTA_CALCULATE_DP(COORDINATE_SYSTEM,X,Y,ERR,ERROR)
!Argument variables
TYPE(COORDINATE_SYSTEM_TYPE), INTENT(IN) :: COORDINATE_SYSTEM !<The coordinate system to calculate the delta for
REAL(DP), INTENT(IN) :: X(:) !<The first coordinate
REAL(DP), INTENT(IN) :: Y(:) !<The second coordinate
INTEGER(INTG), INTENT(OUT) :: ERR !<The error code
TYPE(VARYING_STRING), INTENT(OUT) :: ERROR !<The error string
!Function variable
REAL(DP) :: COORDINATE_DELTA_CALCULATE_DP(SIZE(X,1))
!Local variables
ENTERS("COORDINATE_DELTA_CALCULATE_DP",ERR,ERROR,*999)
COORDINATE_DELTA_CALCULATE_DP=0.0_DP
IF(SIZE(X,1)<COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS) &
& CALL FlagError("Size of X is less than the number of dimensions.",ERR,ERROR,*999)
IF(SIZE(X,1)/=SIZE(Y,1)) &
& CALL FlagError("Size of X is different to the size of Y.",ERR,ERROR,*999)
COORDINATE_DELTA_CALCULATE_DP(1:COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS)=Y(1:COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS)- &
& X(1:COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS)
SELECT CASE(COORDINATE_SYSTEM%TYPE)
CASE(COORDINATE_RECTANGULAR_CARTESIAN_TYPE)
!Do nothing
CASE(COORDINATE_CYLINDRICAL_POLAR_TYPE)
CALL FlagError("Not implemented.",ERR,ERROR,*999)
CASE(COORDINATE_SPHERICAL_POLAR_TYPE)
CALL FlagError("Not implemented.",ERR,ERROR,*999)
CASE(COORDINATE_PROLATE_SPHEROIDAL_TYPE)
CALL FlagError("Not implemented.",ERR,ERROR,*999)
CASE(COORDINATE_OBLATE_SPHEROIDAL_TYPE)
CALL FlagError("Not implemented.",ERR,ERROR,*999)
CASE DEFAULT
CALL FlagError("Invalid coordinate type.",ERR,ERROR,*999)
END SELECT
EXITS("COORDINATE_DELTA_CALCULATE_DP")
RETURN
999 ERRORSEXITS("COORDINATE_DELTA_CALCULATE_DP",ERR,ERROR)
RETURN
END FUNCTION COORDINATE_DELTA_CALCULATE_DP
!
!================================================================================================================================
!
!>Calculates the covariant metric tensor GL(i,j), the contravariant metric tensor GU(i,J), the Jacobian and derivative of the interpolated coordinate system (XI_i) with respect to the given coordinate (X_j) system (DXI_DX) at a point (X - normally a Gauss point). Old cmiss name: XGMG
SUBROUTINE COORDINATE_METRICS_CALCULATE(COORDINATE_SYSTEM,JACOBIAN_TYPE,METRICS,ERR,ERROR,*)
!Argument variables
TYPE(COORDINATE_SYSTEM_TYPE), POINTER :: COORDINATE_SYSTEM !<A pointer to the coordinate system to calculate the metrics for
INTEGER(INTG), INTENT(IN) :: JACOBIAN_TYPE !<The type of Jacobian to calculate \see COORDINATE_ROUTINES_JacobianTypes,COORDINATE_ROUTINES
TYPE(FIELD_INTERPOLATED_POINT_METRICS_TYPE), POINTER :: METRICS !<A pointer to the metrics to calculate
INTEGER(INTG), INTENT(OUT) :: ERR !<The error code
TYPE(VARYING_STRING), INTENT(OUT) :: ERROR !<The error string
!Local Variables
INTEGER(INTG) :: mi,ni,nu
REAL(DP) :: DET_GL,DET_DX_DXI,DX_DXI2(3),DX_DXI3(3),FF,G1,G3,LENGTH,MU,R,RC,RCRC,RR,SCALE
TYPE(FIELD_INTERPOLATED_POINT_TYPE), POINTER :: INTERPOLATED_POINT
TYPE(VARYING_STRING) :: LOCAL_ERROR
NULLIFY(INTERPOLATED_POINT)
ENTERS("COORDINATE_METRICS_CALCULATE",ERR,ERROR,*999)
IF(ASSOCIATED(COORDINATE_SYSTEM)) THEN
IF(ASSOCIATED(METRICS)) THEN
INTERPOLATED_POINT=>METRICS%INTERPOLATED_POINT
IF(ASSOCIATED(INTERPOLATED_POINT)) THEN
IF(INTERPOLATED_POINT%PARTIAL_DERIVATIVE_TYPE>=FIRST_PART_DERIV) THEN
SELECT CASE(METRICS%NUMBER_OF_XI_DIMENSIONS)
CASE(1)
!Calculate the derivatives of X with respect to XI
nu=PARTIAL_DERIVATIVE_FIRST_DERIVATIVE_MAP(1)
METRICS%DX_DXI(1:METRICS%NUMBER_OF_X_DIMENSIONS,1)=INTERPOLATED_POINT%VALUES(1:METRICS%NUMBER_OF_X_DIMENSIONS,nu)
!Initialise the covariant metric tensor to the identity matrix
METRICS%GL(1,1)=1.0_DP
CASE(2)
!Calculate the derivatives of X with respect to XI
nu=PARTIAL_DERIVATIVE_FIRST_DERIVATIVE_MAP(1)
METRICS%DX_DXI(1:METRICS%NUMBER_OF_X_DIMENSIONS,1)=INTERPOLATED_POINT%VALUES(1:METRICS%NUMBER_OF_X_DIMENSIONS,nu)
nu=PARTIAL_DERIVATIVE_FIRST_DERIVATIVE_MAP(2)
METRICS%DX_DXI(1:METRICS%NUMBER_OF_X_DIMENSIONS,2)=INTERPOLATED_POINT%VALUES(1:METRICS%NUMBER_OF_X_DIMENSIONS,nu)
!Initialise the covariant metric tensor to the identity matrix
METRICS%GL(1,1)=1.0_DP
METRICS%GL(1,2)=0.0_DP
METRICS%GL(2,1)=0.0_DP
METRICS%GL(2,2)=1.0_DP
CASE(3)
!Calculate the derivatives of X with respect to XI
nu=PARTIAL_DERIVATIVE_FIRST_DERIVATIVE_MAP(1)
METRICS%DX_DXI(1:METRICS%NUMBER_OF_X_DIMENSIONS,1)=INTERPOLATED_POINT%VALUES(1:METRICS%NUMBER_OF_X_DIMENSIONS,nu)
nu=PARTIAL_DERIVATIVE_FIRST_DERIVATIVE_MAP(2)
METRICS%DX_DXI(1:METRICS%NUMBER_OF_X_DIMENSIONS,2)=INTERPOLATED_POINT%VALUES(1:METRICS%NUMBER_OF_X_DIMENSIONS,nu)
nu=PARTIAL_DERIVATIVE_FIRST_DERIVATIVE_MAP(3)
METRICS%DX_DXI(1:METRICS%NUMBER_OF_X_DIMENSIONS,3)=INTERPOLATED_POINT%VALUES(1:METRICS%NUMBER_OF_X_DIMENSIONS,nu)
!Initialise the covariant metric tensor to the identity matrix
METRICS%GL(1,1)=1.0_DP
METRICS%GL(1,2)=0.0_DP
METRICS%GL(1,3)=0.0_DP
METRICS%GL(2,1)=0.0_DP
METRICS%GL(2,2)=1.0_DP
METRICS%GL(2,3)=0.0_DP
METRICS%GL(3,1)=0.0_DP
METRICS%GL(3,2)=0.0_DP
METRICS%GL(3,3)=1.0_DP
CASE DEFAULT
CALL FlagError("Not implemented.",ERR,ERROR,*999)
END SELECT
!Calculate the covariant metric tensor GL(i,j)
SELECT CASE(COORDINATE_SYSTEM%TYPE)
CASE(COORDINATE_RECTANGULAR_CARTESIAN_TYPE)
SELECT CASE(METRICS%NUMBER_OF_X_DIMENSIONS)
CASE(1)
DO mi=1,METRICS%NUMBER_OF_XI_DIMENSIONS
DO ni=1,METRICS%NUMBER_OF_XI_DIMENSIONS
METRICS%GL(mi,ni)=METRICS%DX_DXI(1,mi)*METRICS%DX_DXI(1,ni)
ENDDO !ni
ENDDO !mi
CASE(2)
DO mi=1,METRICS%NUMBER_OF_XI_DIMENSIONS
DO ni=1,METRICS%NUMBER_OF_XI_DIMENSIONS
METRICS%GL(mi,ni)=METRICS%DX_DXI(1,mi)*METRICS%DX_DXI(1,ni)+ &
& METRICS%DX_DXI(2,mi)*METRICS%DX_DXI(2,ni)
ENDDO !ni
ENDDO !mi
CASE(3)
DO mi=1,METRICS%NUMBER_OF_XI_DIMENSIONS
DO ni=1,METRICS%NUMBER_OF_XI_DIMENSIONS
METRICS%GL(mi,ni)=METRICS%DX_DXI(1,mi)*METRICS%DX_DXI(1,ni)+ &
& METRICS%DX_DXI(2,mi)*METRICS%DX_DXI(2,ni)+ &
& METRICS%DX_DXI(3,mi)*METRICS%DX_DXI(3,ni)
ENDDO !ni
ENDDO !mi
CASE DEFAULT
LOCAL_ERROR=TRIM(NUMBER_TO_VSTRING(METRICS%NUMBER_OF_X_DIMENSIONS,"*",ERR,ERROR))// &
& " is an invalid number of dimensions for a rectangular cartesian coordinate system."
CALL FlagError(LOCAL_ERROR,ERR,ERROR,*999)
END SELECT
CASE(COORDINATE_CYLINDRICAL_POLAR_TYPE)
R=INTERPOLATED_POINT%VALUES(1,1)
RR=R*R
IF(METRICS%NUMBER_OF_X_DIMENSIONS==2) THEN
DO mi=1,METRICS%NUMBER_OF_XI_DIMENSIONS
DO ni=1,METRICS%NUMBER_OF_XI_DIMENSIONS
METRICS%GL(mi,ni)=METRICS%DX_DXI(1,mi)*METRICS%DX_DXI(1,ni)+RR*METRICS%DX_DXI(2,mi)*METRICS%DX_DXI(2,ni)
ENDDO !ni
ENDDO !mi
ELSE IF(METRICS%NUMBER_OF_X_DIMENSIONS==3) THEN
DO mi=1,METRICS%NUMBER_OF_XI_DIMENSIONS
DO ni=1,METRICS%NUMBER_OF_XI_DIMENSIONS
METRICS%GL(mi,ni)=METRICS%DX_DXI(1,mi)*METRICS%DX_DXI(1,ni)+RR*METRICS%DX_DXI(2,mi)*METRICS%DX_DXI(2,ni)+ &
& METRICS%DX_DXI(3,mi)*METRICS%DX_DXI(3,ni)
ENDDO !ni
ENDDO !mi
ELSE
LOCAL_ERROR=TRIM(NUMBER_TO_VSTRING(METRICS%NUMBER_OF_X_DIMENSIONS,"*",ERR,ERROR))// &
& " is an invalid number of dimensions for a cylindrical polar coordinate system."
CALL FlagError(LOCAL_ERROR,ERR,ERROR,*999)
ENDIF
CASE(COORDINATE_SPHERICAL_POLAR_TYPE)
R=INTERPOLATED_POINT%VALUES(1,1)
RR=R*R
RC=R*COS(INTERPOLATED_POINT%VALUES(3,1))
RCRC=RC*RC
DO mi=1,METRICS%NUMBER_OF_XI_DIMENSIONS
DO ni=1,METRICS%NUMBER_OF_XI_DIMENSIONS
METRICS%GL(mi,ni)=METRICS%DX_DXI(1,mi)*METRICS%DX_DXI(1,ni)+RCRC*METRICS%DX_DXI(2,mi)*METRICS%DX_DXI(2,ni)+ &
& RR*METRICS%DX_DXI(3,mi)*METRICS%DX_DXI(3,ni)
ENDDO !ni
ENDDO !mi
CASE(COORDINATE_PROLATE_SPHEROIDAL_TYPE)
IF(ABS(INTERPOLATED_POINT%VALUES(2,1))<ZERO_TOLERANCE) THEN
CALL FLAG_WARNING("Mu is zero.",ERR,ERROR,*999)
ELSE
FF=COORDINATE_SYSTEM%FOCUS*COORDINATE_SYSTEM%FOCUS
R=INTERPOLATED_POINT%VALUES(1,1)
MU=INTERPOLATED_POINT%VALUES(2,1)
G1=FF*(SINH(R)*SINH(R)+SIN(MU)*SIN(MU))
G3=FF*SINH(R)*SINH(R)*SIN(MU)*SIN(MU)
IF(METRICS%NUMBER_OF_X_DIMENSIONS==2) THEN
DO mi=1,METRICS%NUMBER_OF_XI_DIMENSIONS
DO ni=1,METRICS%NUMBER_OF_XI_DIMENSIONS
METRICS%GL(mi,ni)=G1*(METRICS%DX_DXI(1,mi)*METRICS%DX_DXI(1,ni)+METRICS%DX_DXI(2,mi)*METRICS%DX_DXI(2,ni))
ENDDO !ni
ENDDO !mi
ELSE IF(METRICS%NUMBER_OF_X_DIMENSIONS==3) THEN
DO mi=1,METRICS%NUMBER_OF_XI_DIMENSIONS
DO ni=1,METRICS%NUMBER_OF_XI_DIMENSIONS
METRICS%GL(mi,ni)=G1*(METRICS%DX_DXI(1,mi)*METRICS%DX_DXI(1,ni)+METRICS%DX_DXI(2,mi)*METRICS%DX_DXI(2,ni))+ &
& G3*METRICS%DX_DXI(3,mi)*METRICS%DX_DXI(3,ni)
ENDDO !ni
ENDDO !mi
ELSE
LOCAL_ERROR=TRIM(NUMBER_TO_VSTRING(METRICS%NUMBER_OF_X_DIMENSIONS,"*",ERR,ERROR))// &
& " is an invalid number of dimensions for a prolate spheroidal coordinate system."
CALL FlagError(LOCAL_ERROR,ERR,ERROR,*999)
ENDIF
ENDIF
CASE(COORDINATE_OBLATE_SPHEROIDAL_TYPE)
CALL FlagError("Not implemented.",ERR,ERROR,*999)
CASE DEFAULT
LOCAL_ERROR="The coordinate system type of "//TRIM(NUMBER_TO_VSTRING(COORDINATE_SYSTEM%TYPE,"*",ERR,ERROR))// &
& " is invalid."
CALL FlagError(LOCAL_ERROR,ERR,ERROR,*999)
END SELECT
!Calcualte the contravariant metric tensor
CALL INVERT(METRICS%GL(1:METRICS%NUMBER_OF_XI_DIMENSIONS,1:METRICS%NUMBER_OF_XI_DIMENSIONS), &
& METRICS%GU(1:METRICS%NUMBER_OF_XI_DIMENSIONS,1:METRICS%NUMBER_OF_XI_DIMENSIONS),DET_GL, &
& ERR,ERROR,*999)
!Calculate the Jacobian
SELECT CASE(JACOBIAN_TYPE)
CASE(COORDINATE_JACOBIAN_NO_TYPE)
METRICS%JACOBIAN=0.0
METRICS%JACOBIAN_TYPE=COORDINATE_JACOBIAN_NO_TYPE
CASE(COORDINATE_JACOBIAN_LINE_TYPE)
METRICS%JACOBIAN=SQRT(ABS(METRICS%GL(1,1)))
METRICS%JACOBIAN_TYPE=COORDINATE_JACOBIAN_LINE_TYPE
CASE(COORDINATE_JACOBIAN_AREA_TYPE)
IF(METRICS%NUMBER_OF_XI_DIMENSIONS==3) THEN
METRICS%JACOBIAN=SQRT(ABS(DET_GL*METRICS%GU(3,3)))
ELSE
METRICS%JACOBIAN=SQRT(ABS(DET_GL))
ENDIF
METRICS%JACOBIAN_TYPE=COORDINATE_JACOBIAN_AREA_TYPE
CASE(COORDINATE_JACOBIAN_VOLUME_TYPE)
METRICS%JACOBIAN=SQRT(ABS(DET_GL))
METRICS%JACOBIAN_TYPE=COORDINATE_JACOBIAN_VOLUME_TYPE
CASE DEFAULT
LOCAL_ERROR="The Jacobian type of "//TRIM(NUMBER_TO_VSTRING(JACOBIAN_TYPE,"*",ERR,ERROR))// &
& " is invalid."
CALL FlagError(LOCAL_ERROR,ERR,ERROR,*999)
END SELECT
!Calculate the derivatives of Xi with respect to X - DXI_DX
IF(METRICS%NUMBER_OF_XI_DIMENSIONS==METRICS%NUMBER_OF_X_DIMENSIONS) THEN
CALL INVERT(METRICS%DX_DXI,METRICS%DXI_DX,DET_DX_DXI,ERR,ERROR,*999)
ELSE
!We have a line or a surface embedded in a higher dimensional space
SELECT CASE(METRICS%NUMBER_OF_XI_DIMENSIONS)
CASE(1)
!Line in space
SELECT CASE(METRICS%NUMBER_OF_X_DIMENSIONS)
CASE(2)
IF(INTERPOLATED_POINT%PARTIAL_DERIVATIVE_TYPE>FIRST_PART_DERIV) THEN
!We have curvature information. Form the frenet vector frame.
!Calculate the normal vector from the normalised second derivative of the position vector.
nu=PARTIAL_DERIVATIVE_SECOND_DERIVATIVE_MAP(1)
CALL Normalise(INTERPOLATED_POINT%VALUES(1:2,nu),DX_DXI2,ERR,ERROR,*999)
ELSE
!No curvature information but obtain other normal frenet vector by rotating tangent vector 90 deg.
DX_DXI2(1)=-1.0_DP*METRICS%DX_DXI(2,1)
DX_DXI2(2)=METRICS%DX_DXI(1,1)
ENDIF
DET_DX_DXI=METRICS%DX_DXI(1,1)*DX_DXI2(2)-METRICS%DX_DXI(2,1)*DX_DXI2(1)
IF(ABS(DET_DX_DXI)>ZERO_TOLERANCE) THEN
METRICS%DXI_DX(1,1)=DX_DXI2(2)/DET_DX_DXI
METRICS%DXI_DX(1,2)=-1.0_DP*DX_DXI2(1)/DET_DX_DXI
!Normalise to ensure that g^11=g^1.g^1
CALL L2Norm(METRICS%DXI_DX(1,1:2),LENGTH,err,error,*999)
SCALE=SQRT(ABS(METRICS%GU(1,1)))/LENGTH
METRICS%DXI_DX(1,1:2)=SCALE*METRICS%DXI_DX(1,1:2)
ELSE
CALL FLAG_WARNING("Zero determinant. Unable to obtain dxi/dx.",ERR,ERROR,*999)
METRICS%DXI_DX=0.0_DP
ENDIF
CASE(3)
IF(INTERPOLATED_POINT%PARTIAL_DERIVATIVE_TYPE>FIRST_PART_DERIV) THEN
!We have curvature information. Form the frenet vector frame.
!Calculate the normal vector from the normalised second derivative of the position vector.
nu=PARTIAL_DERIVATIVE_SECOND_DERIVATIVE_MAP(1)
CALL Normalise(INTERPOLATED_POINT%VALUES(1:3,nu),DX_DXI2,ERR,ERROR,*999)
!Calculate the bi-normal vector from the normalised cross product of the tangent and normal vectors
CALL NormaliseCrossProduct(METRICS%DX_DXI(1:3,1),DX_DXI2,DX_DXI3,ERR,ERROR,*999)
DET_DX_DXI=METRICS%DX_DXI(1,1)*(DX_DXI2(2)*DX_DXI3(3)-DX_DXI2(3)*DX_DXI3(2))+ &
& DX_DXI2(1)*(METRICS%DX_DXI(3,1)*DX_DXI3(2)-DX_DXI3(3)*METRICS%DX_DXI(2,1))+ &
& DX_DXI3(1)*(METRICS%DX_DXI(2,1)*DX_DXI2(3)-METRICS%DX_DXI(3,1)*DX_DXI2(2))
IF(ABS(DET_DX_DXI)>ZERO_TOLERANCE) THEN
METRICS%DXI_DX(1,1)=(DX_DXI3(3)*DX_DXI2(2)-DX_DXI2(3)*DX_DXI3(2))/DET_DX_DXI
METRICS%DXI_DX(1,2)=-1.0_DP*(DX_DXI3(3)*DX_DXI2(1)-DX_DXI2(3)*DX_DXI3(1))/DET_DX_DXI
METRICS%DXI_DX(1,3)=(DX_DXI3(2)*DX_DXI2(1)-DX_DXI2(2)*DX_DXI3(1))/DET_DX_DXI
!Normalise to ensure that g^11=g^1.g^1
CALL L2Norm(METRICS%DXI_DX(1,1:3),LENGTH,err,error,*999)
SCALE=SQRT(ABS(METRICS%GU(1,1)))/LENGTH
METRICS%DXI_DX(1,1:3)=SCALE*METRICS%DX_DXI(1,1:3)
ELSE
CALL FLAG_WARNING("Zero determinant. Unable to obtain dxi/dx.",ERR,ERROR,*999)
METRICS%DXI_DX=0.0_DP
ENDIF
ELSE
METRICS%DXI_DX(1,1)=METRICS%DX_DXI(1,1)
METRICS%DXI_DX(1,2)=METRICS%DX_DXI(2,1)
METRICS%DXI_DX(1,3)=METRICS%DX_DXI(3,1)
!Normalise to ensure that g^11=g^1.g^1
CALL L2Norm(METRICS%DXI_DX(1,1:3),LENGTH,err,error,*999)
SCALE=SQRT(ABS(METRICS%GU(1,1)))/LENGTH
METRICS%DXI_DX(1,1:3)=SCALE*METRICS%DXI_DX(1,1:3)
ENDIF
CASE DEFAULT
CALL FlagError("Invalid embedding of a line in space.",ERR,ERROR,*999)
END SELECT
CASE(2)
!Surface in space
IF(METRICS%NUMBER_OF_X_DIMENSIONS==3) THEN
!Surface in 3D space.
!Calculate the covariant vectors g^1 and g^2. These are calculated as follows:
!First define g_3=g_1 x g_2, and then define g^1=((g_2 x g_3)_b)/DET_GL and g^2=((g_3 x g_1)_b)/DET_GL.
!The _b means lowering the index with the metric tensor of the curvilinear coordinate system.
!This way we have a consistent set of covariant and covariant vectors, i.e. <g_M,g^N>=delta_M^N.
METRICS%DXI_DX(1,1:3)=(METRICS%GL(2,2)*METRICS%DX_DXI(1:3,1)-METRICS%GL(1,2)*METRICS%DX_DXI(1:3,2))/DET_GL
METRICS%DXI_DX(2,1:3)=(METRICS%GL(1,1)*METRICS%DX_DXI(1:3,2)-METRICS%GL(2,1)*METRICS%DX_DXI(1:3,1))/DET_GL
SELECT CASE(COORDINATE_SYSTEM%TYPE)
CASE(COORDINATE_RECTANGULAR_CARTESIAN_TYPE)
!Do nothing
CASE(COORDINATE_CYLINDRICAL_POLAR_TYPE)
R=INTERPOLATED_POINT%VALUES(1,1)
RR=R*R
METRICS%DXI_DX(1:2,2)=METRICS%DXI_DX(1:2,2)*RR
CASE(COORDINATE_SPHERICAL_POLAR_TYPE)
R=INTERPOLATED_POINT%VALUES(1,1)
RR=R*R
RC=R*COS(INTERPOLATED_POINT%VALUES(3,1))
RCRC=RC*RC
METRICS%DXI_DX(1:2,2)=METRICS%DXI_DX(1:2,2)*RCRC
METRICS%DXI_DX(1:2,3)=METRICS%DXI_DX(1:2,3)*RR
CASE(COORDINATE_PROLATE_SPHEROIDAL_TYPE)
IF(ABS(INTERPOLATED_POINT%VALUES(2,1))<ZERO_TOLERANCE) THEN
CALL FLAG_WARNING("Mu is zero.",ERR,ERROR,*999)
ELSE
FF=COORDINATE_SYSTEM%FOCUS*COORDINATE_SYSTEM%FOCUS
R=INTERPOLATED_POINT%VALUES(1,1)
MU=INTERPOLATED_POINT%VALUES(2,1)
G1=FF*(SINH(R)*SINH(R)+SIN(MU)*SIN(MU))
G3=FF*SINH(R)*SINH(R)*SIN(MU)*SIN(MU)
METRICS%DXI_DX(1:2,1)=METRICS%DXI_DX(1:2,1)*G1
METRICS%DXI_DX(1:2,2)=METRICS%DXI_DX(1:2,2)*G1
METRICS%DXI_DX(1:2,3)=METRICS%DXI_DX(1:2,3)*G3
ENDIF
CASE(COORDINATE_OBLATE_SPHEROIDAL_TYPE)
CALL FlagError("Not implemented.",ERR,ERROR,*999)
CASE DEFAULT
LOCAL_ERROR="The coordinate system type of "//TRIM(NUMBER_TO_VSTRING(COORDINATE_SYSTEM%TYPE,"*",ERR,ERROR))// &
& " is invalid."
CALL FlagError(LOCAL_ERROR,ERR,ERROR,*999)
END SELECT
ELSE
CALL FlagError("Invalid embedding of a surface in space.",ERR,ERROR,*999)
ENDIF
CASE DEFAULT
CALL FlagError("Invalid embedding in space.",ERR,ERROR,*999)
END SELECT
ENDIF
ELSE
CALL FlagError("Metrics interpolated point has not been interpolated to include first derivatives.",ERR,ERROR,*999)
ENDIF
ELSE
CALL FlagError("Metrics interpolated point is not associated.",ERR,ERROR,*999)
ENDIF
ELSE
CALL FlagError("Metrics is not associated.",ERR,ERROR,*999)
ENDIF
ELSE
CALL FlagError("Coordinate system is not associated.",ERR,ERROR,*999)
ENDIF
IF(DIAGNOSTICS1) THEN
CALL WriteString(DIAGNOSTIC_OUTPUT_TYPE,"",err,error,*999)
CALL WRITE_STRING(DIAGNOSTIC_OUTPUT_TYPE,"Coordinate system metrics:",ERR,ERROR,*999)
CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE," Coordinate system type = ",TRIM(COORDINATE_SYSTEM_TYPE_STRING( &
& COORDINATE_SYSTEM%TYPE)),ERR,ERROR,*999)
CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE," Number of X dimensions = ",METRICS%NUMBER_OF_X_DIMENSIONS,ERR,ERROR,*999)
CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE," Number of Xi dimensions = ",METRICS%NUMBER_OF_XI_DIMENSIONS,ERR,ERROR,*999)
CALL WRITE_STRING(DIAGNOSTIC_OUTPUT_TYPE," Location of metrics:",ERR,ERROR,*999)
CALL WRITE_STRING_VECTOR(DIAGNOSTIC_OUTPUT_TYPE,1,1,METRICS%NUMBER_OF_X_DIMENSIONS,3,3,INTERPOLATED_POINT%VALUES(:,1), &
& '(" X :",3(X,E13.6))','(17X,3(X,E13.6))',ERR,ERROR,*999)
CALL WRITE_STRING(DIAGNOSTIC_OUTPUT_TYPE," Derivative of X wrt Xi:",ERR,ERROR,*999)
CALL WRITE_STRING_MATRIX(DIAGNOSTIC_OUTPUT_TYPE,1,1,METRICS%NUMBER_OF_X_DIMENSIONS,1,1,METRICS%NUMBER_OF_XI_DIMENSIONS, &
& 3,3,METRICS%DX_DXI,WRITE_STRING_MATRIX_NAME_AND_INDICES,'(" dX_dXi','(",I1,",:)',' :",3(X,E13.6))', &
& '(17X,3(X,E13.6))',ERR,ERROR,*999)
IF(METRICS%NUMBER_OF_X_DIMENSIONS/=METRICS%NUMBER_OF_XI_DIMENSIONS) THEN
CALL WRITE_STRING(DIAGNOSTIC_OUTPUT_TYPE," Constructed derivative of X wrt Xi:",ERR,ERROR,*999)
SELECT CASE(METRICS%NUMBER_OF_XI_DIMENSIONS)
CASE(1)
!Line in space
SELECT CASE(METRICS%NUMBER_OF_X_DIMENSIONS)
CASE(2)
CALL WRITE_STRING_VECTOR(DIAGNOSTIC_OUTPUT_TYPE,1,1,METRICS%NUMBER_OF_X_DIMENSIONS,3,3,DX_DXI2, &
& '(" dX_dXi(:,2) :",3(X,E13.6))','(17X,3(X,E13.6))',ERR,ERROR,*999)
CASE(3)
CALL WRITE_STRING_VECTOR(DIAGNOSTIC_OUTPUT_TYPE,1,1,METRICS%NUMBER_OF_X_DIMENSIONS,3,3,DX_DXI2, &
& '(" dX_dXi(:,2) :",3(X,E13.6))','(17X,3(X,E13.6))',ERR,ERROR,*999)
CALL WRITE_STRING_VECTOR(DIAGNOSTIC_OUTPUT_TYPE,1,1,METRICS%NUMBER_OF_X_DIMENSIONS,3,3,DX_DXI3, &
& '(" dX_dXi(:,3) :",3(X,E13.6))','(17X,3(X,E13.6))',ERR,ERROR,*999)
CASE DEFAULT
CALL FlagError("Invalid embedding of a line in space.",ERR,ERROR,*999)
END SELECT
CASE(2)
!Surface in space
SELECT CASE(METRICS%NUMBER_OF_X_DIMENSIONS)
CASE(3)
CALL WRITE_STRING_VECTOR(DIAGNOSTIC_OUTPUT_TYPE,1,1,METRICS%NUMBER_OF_X_DIMENSIONS,3,3,DX_DXI3, &
& '(" dX_dXi(:,3) :",3(X,E13.6))','(17X,3(X,E13.6))',ERR,ERROR,*999)
CASE DEFAULT
CALL FlagError("Invalid embedding of a surface in space.",ERR,ERROR,*999)
END SELECT