Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Resolve y-axis and w based local coordinate system #440

Merged
merged 1 commit into from
Mar 15, 2024
Merged

Resolve y-axis and w based local coordinate system #440

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
99 changes: 66 additions & 33 deletions fem/src/modules/CoordinateTransform.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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]

Expand Down Expand Up @@ -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]
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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

Expand All @@ -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

Expand Down
74 changes: 74 additions & 0 deletions fem/tests/circuits_harmonic_homogenization_coil_solver/2241/entities.sif
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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

Loading