Resolve y-axis and w based local coordinate system

The main application for the y-axis and w based local coordinate system
is the stranded homogenization model where the homogenization parameters
are oriented using the RotM transformation. The RotM transformation
needs a local coordinate system and the wish is to be able to define it
without direction solvers (in case of arbitrary wire cross-sections).

A test case derived from `circuits_harmonic_stranded_homogenization` is
added. The original case computes the skin and proximity losses for
a stranded coilcurrent via so called homogenization, where the strands
are not geometrically in the model but their frequency responce is
estimated using `sigma` and `nu` parameters that are gotten from the
elementary solutions of 2D computations (transversal fields and
perpendicular current). In 3D the cross-section needs to be oriented
using the RotM transformation that requires the definition of a local
coordinate system. For that we have previously used the direction
solvers (`Alpha` and `Beta`).

The main point now is to get rid of the direction solvers.
The other major change is to use `CoilSolver` solution to drive
the circuit (instead of `WSolver`). There are many benefits for
using the `CoilSolver`, for example, one can calculate closed coils.

The important new keywords (with respect to ) here are:
 1. `Body: Local Coordinate System Beta Reference and Gamma = Logical True`
   * CoordinateTransform computes alpha = beta x gamma:
     where beta is set to "Beta Reference" and
     Gamma is set to "coilcurrent e" computed by `CoilSolver`
 2. `Component: Coil Normal(3) = Real 0. 0. 1.`
   * `CoilSolver` needs to know how the coil is oriented
 3. `Component: Desired Current Density = Real 1.0`
   * `CoilSolver` will set a certain current density over the coil
 4. `Component: Coil Use W Vector = Logical True`
   * Coil uses a defined `W Vector` for which a name can be provided
   * The `W Vector` set here is used for driving the component with circuits
 5. `Component: W Vector Variable Name = String "CoilCurrent e"`
   * Here we choose the correct field for `W Vector`
 6. `Boundary: Coil Start = Logical True`
   * `CoilSolver` needs a boundary for knowing where coil starts (if not closed)
 7. `Boundary: Coil End = Logical True`
   * `CoilSolver` needs a boundary for knowing where coil ends (if not closed)

Issue #429

fem/src/modules/CoordinateTransform.F90

@@ -172,6 +172,8 @@ SUBROUTINE RotMSolver( Model,Solver,dt,TransientSimulation )
   REAL :: PDDetTol
   INTEGER :: PDMaxIter
 
+  LOGICAL :: LocalSystemBetaRefAndGamma
+
   INTEGER, PARAMETER :: ind1(9) = [1,1,1,2,2,2,3,3,3]
   INTEGER, PARAMETER :: ind2(9) = [1,2,3,1,2,3,1,2,3]
 
@@ -224,7 +226,7 @@ SUBROUTINE RotMSolver( Model,Solver,dt,TransientSimulation )
       CALL Warn('CoordinateTransform','Polar Decomposition Max Iteration not set.') 
       PDMaxIter= 100
     END IF  
-    
+
     CoordSys_ijk(1,:) = [1,0,0]
     CoordSys_ijk(2,:) = [0,1,0]
     CoordSys_ijk(3,:) = [0,0,1]
@@ -249,6 +251,9 @@ SUBROUTINE RotMSolver( Model,Solver,dt,TransientSimulation )
     IF (SIZE(beta_ref,1) /= 3) CALL Fatal('RotMSolver','Beta Reference should have three components!')
     CoordSys_ref(2,1:3) = normalized(beta_ref(1:3,1))
 
+    LocalSystemBetaRefAndGamma = ListGetLogical(BodyParams, 'Local Coordinate System Beta Reference and Gamma', Found)
+    IF (.NOT. Found) LocalSystemBetaRefAndGamma = .FALSE.
+
     CoordSys_ref(3,1:3) = normalized(crossproduct(CoordSys_ref(1,1:3), CoordSys_ref(2,1:3)))
 
     ! Compute the rotation matrices for 2 rank tensors for all the elements
@@ -270,7 +275,7 @@ SUBROUTINE ComputeRotM(Element,CoordSys_ijk,CoordSys_ref,nn,nd, &
     REAL(KIND=dp) :: un,vn,wn,Basis(nn), DetJ,localC,gradv(3)
     REAL(KIND=dp) :: dBasisdx(nn,3),Coordsys(3,3),Coordsys2(3,3), &
                      Relem(3,3), jac_ref(3,3), jac_e(3,3), alpha(nn), &
-                     beta(nn), tvec(3)
+                     beta(nn), tvec(3), gamma(3,nn)
     REAL(KIND=dp) :: CoordSys_ijk(3,3), CoordSys_ref(3,3)
 
     INTEGER :: j,Indexes(nd),l,m,k
@@ -282,23 +287,34 @@ SUBROUTINE ComputeRotM(Element,CoordSys_ijk,CoordSys_ref,nn,nd, &
     INTEGER :: PDMaxIter
 
     CALL GetElementNodes(Nodes)
-    tmpvar => VariableGet( Mesh % Variables, 'alpha direction')
-    IF(ASSOCIATED(tmpvar)) THEN
-      CALL GetLocalSolution(alpha,'alpha direction')
+
+    IF (LocalSystemBetaRefAndGamma) THEN
+      gamma=0._dp
+      tmpvar => VariableGet( Mesh % Variables, 'coilcurrent e')
+      IF(ASSOCIATED(tmpvar)) THEN
+        CALL GetLocalSolution(gamma,'coilcurrent e')
+      ELSE
+        CALL Fatal('ComputeRotM', 'Local System Beta Reference And Gamma tried but Gamma (coilcurrent e) is not found!')
+      END IF
     ELSE
-      tmpvar => VariableGet( Mesh % Variables, 'alpha')
+      tmpvar => VariableGet( Mesh % Variables, 'alpha direction')
       IF(ASSOCIATED(tmpvar)) THEN
-        CALL GetLocalSolution(alpha,'alpha')
+        CALL GetLocalSolution(alpha,'alpha direction')
+      ELSE
+        tmpvar => VariableGet( Mesh % Variables, 'alpha')
+        IF(ASSOCIATED(tmpvar)) THEN
+          CALL GetLocalSolution(alpha,'alpha')
+        END IF
       END IF
-    END IF
 
-    tmpvar => VariableGet( Mesh % Variables, 'beta direction')
-    IF(ASSOCIATED(tmpvar)) THEN
-      CALL GetLocalSolution(beta,'beta direction')
-    ELSE
-      tmpvar => VariableGet( Mesh % Variables, 'beta')
+      tmpvar => VariableGet( Mesh % Variables, 'beta direction')
       IF(ASSOCIATED(tmpvar)) THEN
-        CALL GetLocalSolution(beta,'beta')
+        CALL GetLocalSolution(beta,'beta direction')
+      ELSE
+        tmpvar => VariableGet( Mesh % Variables, 'beta')
+        IF(ASSOCIATED(tmpvar)) THEN
+          CALL GetLocalSolution(beta,'beta')
+        END IF
       END IF
     END IF
 
@@ -314,27 +330,44 @@ SUBROUTINE ComputeRotM(Element,CoordSys_ijk,CoordSys_ref,nn,nd, &
       ! CoordSys(1,1:3) is the perpendicular direction to the foil 
       ! surface
       ! -----------------------------------------------------
-      CoordSys(1,1:3) = normalized(MATMUL( alpha(1:nn), dBasisdx(1:nn,:)))
-      IF (ANY(ISNAN(CoordSys(1,:)))) THEN
-        print *, "Element index = ", GetElementIndex(Element)
-        print *, "Element aspect ratio = ", ElementAspectRatio(Model, Element)
-        CALL Warn('CoordinateTransform','Element coordinate system is NaN, this could be &
-          due to a poor mesh. Let us try to use the degenerate element normal as the local coordinate system alpha vector.') 
-        CoordSys(1,1:3) = NormalOfDegenerateElement(Model, Element)
-        IF (ANY(ISNAN(CoordSys(1,:)))) CALL Fatal('CoordinateTransform','Degenerate element normal did not work...') 
-      END IF
+      IF (LocalSystemBetaRefAndGamma) THEN
+        CoordSys(3,1:3) = MATMUL(gamma(1:3,1:nn), basis(1:nn)) ! Assume this is from normalized coil current
+        IF (ANY(ISNAN(CoordSys(3,:)))) THEN
+          print *, "Element index = ", GetElementIndex(Element)
+          print *, "Element aspect ratio = ", ElementAspectRatio(Model, Element)
+          CALL Warn('CoordinateTransform','Element coordinate system is NaN, this could be &
+            due to a poor mesh. Let us try to use the degenerate element normal as the local coordinate system alpha vector.') 
+          CoordSys(3,1:3) = NormalOfDegenerateElement(Model, Element)
+          IF (ANY(ISNAN(CoordSys(3,:)))) CALL Fatal('CoordinateTransform','Degenerate element normal did not work...') 
+        END IF
 
-      CoordSys(2,1:3) = normalized(MATMUL( beta(1:nn), dBasisdx(1:nn,:)))
-      IF (ANY(ISNAN(CoordSys(2,:)))) THEN
-        print *, "Element index = ", GetElementIndex(Element)
-        print *, "Element aspect ratio = ", ElementAspectRatio(Model, Element)
-        CALL Warn('CoordinateTransform','Element coordinate system is NaN, this could be &
-          due to a poor mesh. Let us try to use the degenerate element normal as the local coordinate system beta vector.') 
-        CoordSys(2,1:3) = NormalOfDegenerateElement(Model, Element)
-        IF (ANY(ISNAN(CoordSys(2,:)))) CALL Fatal('CoordinateTransform','Degenerate element normal did not work...') 
-      END IF
+        CoordSys(2,1:3) = CoordSys_ref(2,1:3)
 
-      CoordSys(3,1:3) = normalized(crossproduct(CoordSys(1,1:3), CoordSys(2,1:3)))
+        CoordSys(1,1:3) = crossproduct(CoordSys(2,1:3), CoordSys(3,1:3))
+!        print "('>',3(F5.1,x),/,x)", CoordSys
+      ELSE
+        CoordSys(1,1:3) = normalized(MATMUL( alpha(1:nn), dBasisdx(1:nn,:)))
+        IF (ANY(ISNAN(CoordSys(1,:)))) THEN
+          print *, "Element index = ", GetElementIndex(Element)
+          print *, "Element aspect ratio = ", ElementAspectRatio(Model, Element)
+          CALL Warn('CoordinateTransform','Element coordinate system is NaN, this could be &
+            due to a poor mesh. Let us try to use the degenerate element normal as the local coordinate system alpha vector.') 
+          CoordSys(1,1:3) = NormalOfDegenerateElement(Model, Element)
+          IF (ANY(ISNAN(CoordSys(1,:)))) CALL Fatal('CoordinateTransform','Degenerate element normal did not work...') 
+        END IF
+
+        CoordSys(2,1:3) = normalized(MATMUL( beta(1:nn), dBasisdx(1:nn,:)))
+        IF (ANY(ISNAN(CoordSys(2,:)))) THEN
+          print *, "Element index = ", GetElementIndex(Element)
+          print *, "Element aspect ratio = ", ElementAspectRatio(Model, Element)
+          CALL Warn('CoordinateTransform','Element coordinate system is NaN, this could be &
+            due to a poor mesh. Let us try to use the degenerate element normal as the local coordinate system beta vector.') 
+          CoordSys(2,1:3) = NormalOfDegenerateElement(Model, Element)
+          IF (ANY(ISNAN(CoordSys(2,:)))) CALL Fatal('CoordinateTransform','Degenerate element normal did not work...') 
+        END IF
+
+        CoordSys(3,1:3) = normalized(crossproduct(CoordSys(1,1:3), CoordSys(2,1:3)))
+      END IF
 
       CoordSys2 = CoordSys
 

fem/tests/circuits_harmonic_homogenization_coil_solver/2241/entities.sif

@@ -0,0 +1,74 @@
+!------ Skeleton for body section -----
+Body 1
+  Name = core
+End
+
+Body 2
+  Name = airgap_1_limb_2
+End
+
+Body 3
+  Name = wp1
+End
+
+Body 4
+  Name = air
+End
+
+!------ Skeleton for boundary section -----
+Boundary Condition 1
+  Name = xy0
+End
+
+Boundary Condition 2
+  Name = airgap_1_limb_2_boundary0
+End
+
+Boundary Condition 3
+  Name = wp1_alpha0
+End
+
+Boundary Condition 4
+  Name = wp1_beta1
+End
+
+Boundary Condition 5
+  Name = wp1_alpha1
+End
+
+Boundary Condition 6
+  Name = wp1_beta0
+End
+
+Boundary Condition 7
+  Name = airgap_1_limb_2_boundary1
+End
+
+Boundary Condition 8
+  Name = coreface_xy
+End
+
+Boundary Condition 9
+  Name = wp1_gamma0
+End
+
+Boundary Condition 10
+  Name = cylinder_lateral_surface
+End
+
+Boundary Condition 11
+  Name = wp1_gamma1
+End
+
+Boundary Condition 12
+  Name = infinity_xz1
+End
+
+Boundary Condition 13
+  Name = infinity_xz0
+End
+
+Boundary Condition 14
+  Name = bnry18
+End
+