Skip to content

Commit 4bc52cd

Browse files
rabackraback
committed
Critical fix.
1 parent 35fdf63 commit 4bc52cd

File tree

3 files changed

+268
-87
lines changed

3 files changed

+268
-87
lines changed

‎fem/src/MortarUtils.F90‎

Lines changed: 75 additions & 66 deletions
Original file line numberDiff line numberDiff line change
@@ -1580,9 +1580,11 @@ END SUBROUTINE TemporalSegmentMortarAssembly
15801580
!----------------------------------------------------------------------------------------
15811581
!> Given a temporal segment "ElementT", calculate mass matrix contributions for projection
15821582
!> for the slave element "Element" and master element "ElementM".
1583+
!> This routine implements the Nitsche method for the interface condition. Note that
1584+
!> the numbering here is that of the initial mesh, not of the temporal interface mesh!
15831585
!----------------------------------------------------------------------------------------
15841586
SUBROUTINE TemporalSegmentNitscheAssembly(ElementT, NodesT, Element, Nodes, ElementM, NodesM, &
1585-
sgn0, pElemBasis, NoGaussPoints, Projector, NodeCoeff, ArcCoeff, NodeScale, &
1587+
sgn0, pElemBasis, NoGaussPoints, Projector, ArcCoeff, NodeScale, &
15861588
NodePerm, InvPerm, InvPermM, SumArea, BC )
15871589
!----------------------------------------------------------------------------------------
15881590
TYPE(Element_t) :: ElementT
@@ -1592,7 +1594,7 @@ SUBROUTINE TemporalSegmentNitscheAssembly(ElementT, NodesT, Element, Nodes, Elem
15921594
LOGICAL :: pElemBasis
15931595
INTEGER :: NoGaussPoints
15941596
TYPE(Matrix_t) :: Projector
1595-
REAL(KIND=dp) :: NodeCoeff, ArcCoeff, NodeScale, SumArea
1597+
REAL(KIND=dp) :: ArcCoeff, NodeScale, SumArea
15961598
INTEGER :: NodePerm(:)
15971599
INTEGER, ALLOCATABLE :: DualNodePerm(:)
15981600
INTEGER, POINTER :: InvPerm(:), InvPermM(:)
@@ -1603,31 +1605,28 @@ SUBROUTINE TemporalSegmentNitscheAssembly(ElementT, NodesT, Element, Nodes, Elem
16031605
INTEGER, POINTER :: Indexes(:), IndexesM(:)
16041606
REAL(KIND=dp) :: val, val_dual, u, v, w, um, vm, wm, xt, yt, zt, weight,DetJ
16051607
REAL(KIND=dp), ALLOCATABLE :: BasisT(:),Basis(:), BasisM(:),pBasis(:),pBasisM(:),dBasisdx(:,:), dBasisdxM(:,:)
1606-
LOGICAL :: AllocationsDone = .FALSE.
1608+
LOGICAL :: AllocationsDone = .FALSE., Found
16071609
TYPE(Matrix_t), POINTER :: DualProjector
16081610
LOGICAL :: Stat
16091611

16101612
TYPE(Mesh_t), POINTER :: Mesh
16111613
TYPE(Element_t), POINTER :: TrueElement, TrueElementM, Parent, ParentM
16121614
TYPE(Nodes_t), SAVE :: TrueNodes, TrueNodesM, ParentNodes, ParentNodesM
16131615
INTEGER :: np, npM, Lcols(20)
1614-
REAL(KIND=dp) :: Nrm(3), NrmM(3), Esize, EsizeM, Gamma, LVals(20)
1616+
REAL(KIND=dp) :: Nrm(3), NrmM(3), Esize, EsizeM, Gamma, Cond, LVals(20)
16151617
INTEGER, ALLOCATABLE, TARGET :: Ind(:), IndM(:), pIndexes(:), pIndexesM(:)
1616-
INTEGER, SAVE :: sgns(4)
1618+
INTEGER, SAVE :: sgns(4), previ=-1
16171619

16181620

16191621
SAVE :: BasisT, Basis, BasisM, pBasis, pBasisM, dBasisdx, dBasisdxM, Ind, IndM, pIndexes, pIndexesM
16201622

16211623
Mesh => CurrentModel % Mesh
16221624

1623-
IF(.NOT. AllocationsDone ) THEN
1624-
n = CurrentModel % Mesh % MaxElementDofs
1625-
ALLOCATE(BasisT(n),Basis(n),BasisM(n),pBasis(n),pBasisM(n),dBasisdx(n,3), dBasisdxM(n,3), &
1626-
Indexes(n), IndexesM(n), Ind(n), IndM(n))
1627-
1628-
i = ListGetInteger( CurrentModel % Simulation,'signs')
1629-
1630-
sgns = 1
1625+
! Testing to test different sign conventions.
1626+
i = NINT( ListGetCReal( CurrentModel % Simulation,'signs', Found ) )
1627+
sgns = 1
1628+
IF( Found .AND. i /= previ ) THEN
1629+
previ = i
16311630
IF(i>=8) THEN
16321631
sgns(1) = -1
16331632
i=i-8
@@ -1644,9 +1643,14 @@ SUBROUTINE TemporalSegmentNitscheAssembly(ElementT, NodesT, Element, Nodes, Elem
16441643
sgns(4) = -1
16451644
i=i-1
16461645
END IF
1647-
1648-
PRINT *,'signs:',sgns
1649-
1646+
PRINT *,'signs:',previ,sgns
1647+
END IF
1648+
1649+
1650+
IF(.NOT. AllocationsDone ) THEN
1651+
n = CurrentModel % Mesh % MaxElementDofs
1652+
ALLOCATE(BasisT(n),Basis(n),BasisM(n),pBasis(n),pBasisM(n),dBasisdx(n,3), dBasisdxM(n,3), &
1653+
Indexes(n), IndexesM(n), Ind(n), IndM(n))
16501654
AllocationsDone = .TRUE.
16511655
END IF
16521656

@@ -1689,10 +1693,6 @@ SUBROUTINE TemporalSegmentNitscheAssembly(ElementT, NodesT, Element, Nodes, Elem
16891693
Nrm = NormalVector( TrueElement, TrueNodes, Check = .TRUE. )
16901694
NrmM = NormalVector( TrueElementM, TrueNodesM, Check = .TRUE. )
16911695

1692-
! Swap either to get consistency...
1693-
NrmM = sgns(4) * NrmM
1694-
1695-
16961696
Esize = ElementDiameter(TrueElement, TrueNodes)
16971697
EsizeM = ElementDiameter(TrueElementM, TrueNodesM)
16981698

@@ -1702,20 +1702,19 @@ SUBROUTINE TemporalSegmentNitscheAssembly(ElementT, NodesT, Element, Nodes, Elem
17021702
END IF
17031703
ParentM => ElementM % BoundaryInfo % Left
17041704
IF(.NOT. ASSOCIATED(ParentM)) THEN
1705-
ParentM => ElementM% BoundaryInfo % Right
1705+
ParentM => ElementM % BoundaryInfo % Right
17061706
END IF
1707-
1707+
17081708
np = Parent % Type % NumberOfNodes
17091709
npM = ParentM % Type % NumberOfNodes
17101710
CALL CopyElementNodesFromMesh(ParentNodes, Mesh, np, Parent % NodeIndexes)
17111711
CALL CopyElementNodesFromMesh(ParentNodesM, Mesh, npM, ParentM % NodeIndexes)
17121712

1713-
Gamma = ListGetCReal(BC,'Nitsche Penalty')
1713+
Gamma = ListGetCReal(BC,'Nitsche Penalty',Found)
1714+
IF(.NOT. Found) Gamma = 1.0e-3
1715+
Cond = ListGetCReal(BC,'Nitsche Conductivity',Found)
1716+
IF(.NOT. Found) Cond = 1.0_dp
17141717

1715-
1716-
1717-
PRINT *,'Elem:',Nrm, np, Esize
1718-
PRINT *,'ElemM:',NrmM, npm, EsizeM
17191718

17201719
DO nip=1, IPT % n
17211720
stat = ElementInfo( ElementT,NodesT,IPT % u(nip),&
@@ -1730,18 +1729,22 @@ SUBROUTINE TemporalSegmentNitscheAssembly(ElementT, NodesT, Element, Nodes, Elem
17301729
Weight = ArcCoeff * DetJ * IPT % s(nip)
17311730
sumarea = sumarea + Weight
17321731

1732+
! We scale the whole condition with conductivity since all other terms in the stiffness
1733+
! matrix are also scaled with this...
1734+
Weight = Cond * Weight
1735+
17331736
! Slave element
17341737
! Integration point, working on 2D plane
17351738
CALL GlobalToLocal( u, v, w, xt, yt, zt, Element, Nodes )
1736-
stat = ElementInfo( Element, Nodes, u, v, w, detJ, Basis )
1737-
1739+
!stat = ElementInfo( Element, Nodes, u, v, w, detJ, Basis )
1740+
stat = ElementInfo( TrueElement, TrueNodes, u, v, w, detJ, Basis )
1741+
17381742
! Basis functions at parent element
17391743
CALL FindParentUVW( TrueElement, n, Parent, np, U, V, W, Basis )
1740-
stat = ElementInfo( Parent, ParentNodes, &
1741-
U, V, W, detJ, pBasis, dBasisdx )
1744+
stat = ElementInfo( Parent, ParentNodes, U, V, W, detJ, pBasis, dBasisdx )
17421745

17431746
! Mapping from boundary local index to parent local index
1744-
DO i=1,nd
1747+
DO i=1,n
17451748
DO ii=1,np
17461749
IF ( TrueElement % NodeIndexes(i) == Parent % NodeIndexes(ii) ) THEN
17471750
Ind(i) = ii; EXIT
@@ -1751,8 +1754,9 @@ SUBROUTINE TemporalSegmentNitscheAssembly(ElementT, NodesT, Element, Nodes, Elem
17511754

17521755
! Master element, same steps
17531756
CALL GlobalToLocal( um, vm, wm, xt, yt, zt, ElementM, NodesM )
1754-
stat = ElementInfo( ElementM, NodesM, um, vm, wm, detJ, BasisM )
1755-
1757+
!stat = ElementInfo( ElementM, NodesM, um, vm, wm, detJ, BasisM )
1758+
stat = ElementInfo( TrueElementM, TrueNodesM, um, vm, wm, detJ, BasisM )
1759+
17561760
CALL FindParentUVW( TrueElementM, nM, ParentM, npM, Um, Vm, Wm, BasisM )
17571761
stat = ElementInfo( ParentM, ParentNodesM, &
17581762
Um, Vm, Wm, detJ, pBasisM, dBasisdxM )
@@ -1764,64 +1768,62 @@ SUBROUTINE TemporalSegmentNitscheAssembly(ElementT, NodesT, Element, Nodes, Elem
17641768
END IF
17651769
END DO
17661770
END DO
1767-
1771+
17681772
! Add the entries to the projector
17691773
DO j=1,nd
1774+
! Testing with Basis(j)
17701775
jj = Indexes(j)
17711776
IF (j<=n) jj = InvPerm(jj)
17721777

1773-
!Projector % InvPerm(nrow) = jj
1774-
!val = NodeCoeff * Basis(j) * Weight
1775-
17761778
DO i=1,nd
17771779
ii = Indexes(i)
17781780
IF(i<=n) ii = InvPerm(ii)
17791781

17801782
LCols(i) = ii
1781-
LVals(i) = weight * ( SUM(dBasisdx(Ind(j),:)*Nrm) * Basis(i) &
1782-
+ SUM(dBasisdx(Ind(i),:)*Nrm) * Basis(j) &
1783-
+ Basis(i) * Basis(j) / Esize / Gamma )
1783+
LVals(i) = sgns(1) * SUM(dBasisdx(Ind(j),:)*Nrm) * Basis(i) & ! consistency term to weakly enforce continuity of solution
1784+
+ sgns(2) * SUM(dBasisdx(Ind(i),:)*Nrm) * Basis(j) & ! consistency term to weakly enforce continuity of flux
1785+
+ Basis(i) * Basis(j) / Esize / Gamma ! penalty term
17841786
END DO
17851787

17861788
DO i=1,ndM
17871789
ii = IndexesM(i)
1788-
IF(i<=nM) ii=InvPermM(ii)
1790+
IF(i<=nM) ii = InvPermM(ii)
17891791

17901792
LCols(nd+i) = ii
1791-
LVals(nd+i) = -weight * ( sgns(1) * SUM(dBasisdx(Ind(j),:)*Nrm) * BasisM(i) &
1792-
+ sgns(2) * SUM(dBasisdxM(IndM(i),:)*NrmM) * Basis(j) &
1793-
+ sgns(3) * BasisM(i) * Basis(j) / Esize / Gamma )
1793+
LVals(nd+i) = -NodeScale * ( sgns(3) * SUM(dBasisdx(Ind(j),:)*Nrm) * BasisM(i) &
1794+
+ sgns(4) * SUM(dBasisdxM(IndM(i),:)*NrmM) * Basis(j) &
1795+
+ BasisM(i) * Basis(j) / Esize / Gamma )
17941796
END DO
1797+
LVals(1:nd+ndM) = weight * LVals(1:nd+ndM)
17951798
CALL List_AddMatrixRow(Projector % ListMatrix,jj,nd+ndM,Lcols,Lvals)
17961799
END DO
17971800

17981801
! Dual projector
17991802
DO j=1,ndM
18001803
jj = IndexesM(j)
1801-
IF (j<=n) jj = InvPermM(jj)
1802-
1803-
!Projector % InvPerm(nrow) = jj
1804-
!val = NodeCoeff * BasisM(j) * Weight
1804+
IF (j<=nM) jj = InvPermM(jj)
18051805

18061806
DO i=1,ndM
18071807
ii = IndexesM(i)
1808-
IF(i<=n) ii=InvPermM(ii)
1808+
IF(i<=nM) ii = InvPermM(ii)
18091809

18101810
LCols(i) = ii
1811-
LVals(i) = weight * ( SUM(dBasisdxM(Ind(j),:)*NrmM) * BasisM(i) &
1812-
+ SUM(dBasisdxM(Ind(i),:)*NrmM) * BasisM(j) &
1813-
+ BasisM(i) * BasisM(j) / EsizeM / Gamma )
1811+
LVals(i) = sgns(1) * SUM(dBasisdxM(Ind(j),:)*NrmM) * BasisM(i) &
1812+
+ sgns(2) * SUM(dBasisdxM(IndM(i),:)*NrmM) * BasisM(j) &
1813+
+ BasisM(i) * BasisM(j) / EsizeM / Gamma
18141814
END DO
18151815

18161816
DO i=1,nd
18171817
ii = Indexes(i)
1818-
IF(i<=nM) ii=InvPerm(ii)
1818+
IF(i<=n) ii = InvPerm(ii)
18191819

18201820
LCols(ndM+i) = ii
1821-
LVals(ndM+i) = -weight * ( sgns(1) * SUM(dBasisdxM(Ind(j),:)*NrmM) * Basis(i) &
1822-
+ sgns(2) * SUM(dBasisdx(Ind(i),:)*Nrm) * BasisM(j) &
1823-
+ sgns(3) * Basis(i) * BasisM(j) / EsizeM / Gamma )
1824-
END DO
1821+
LVals(ndM+i) = -NodeScale * ( sgns(3) * SUM(dBasisdxM(Ind(j),:)*NrmM) * Basis(i) &
1822+
+ sgns(4) * SUM(dBasisdx(Ind(i),:)*Nrm) * BasisM(j) &
1823+
+ Basis(i) * BasisM(j) / EsizeM / Gamma )
1824+
END DO
1825+
1826+
LVals(1:nd+ndM) = weight * LVals(1:nd+ndM)
18251827
CALL List_AddMatrixRow(Projector % ListMatrix,jj,ndM+nd,Lcols,Lvals)
18261828
END DO
18271829
END DO
@@ -3059,8 +3061,7 @@ SUBROUTINE NormalProjectorWeak1D()
30593061
! restructured into a separate routine.
30603062
IF( Projector % ProjectorType == PROJECTOR_TYPE_NITSCHE ) THEN
30613063
CALL TemporalSegmentNitscheAssembly(ElementT, NodesT, Element, Nodes, ElementM, NodesM, &
3062-
sgn0, pElemBasis, 0, &
3063-
Projector, NodeCoeff, ArcCoeff, NodeScale, NodePerm, &
3064+
sgn0, pElemBasis, 0, Projector, ArcCoeff, NodeScale, NodePerm, &
30643065
InvPerm1, InvPerm2, SumArea, BC )
30653066
ELSE
30663067
CALL TemporalSegmentMortarAssembly(ElementT, NodesT, Element, Nodes, ElementM, NodesM, &
@@ -3309,7 +3310,7 @@ FUNCTION LevelProjector( BMesh1, BMesh2, Repeating, AntiRepeating, &
33093310
REAL(KIND=dp), ALLOCATABLE :: Cond(:)
33103311
TYPE(Matrix_t), POINTER :: DualProjector
33113312
LOGICAL :: DualMaster, DualSlave, DualLCoeff, BiorthogonalBasis
3312-
LOGICAL :: SecondFamily, SecondOrder, pElemProj, pElemBasis
3313+
LOGICAL :: SecondFamily, SecondOrder, pElemProj, pElemBasis, NitscheProjector
33133314
CHARACTER(*), PARAMETER :: Caller = "LevelProjector"
33143315

33153316
CALL Info(Caller,'Creating projector for a levelized mesh',Level=7)
@@ -3378,8 +3379,16 @@ FUNCTION LevelProjector( BMesh1, BMesh2, Repeating, AntiRepeating, &
33783379
ArcCoeff = 1.0_dp
33793380
END IF
33803381

3382+
3383+
33813384
! We have a weak projector if it is requested
3382-
WeakProjector = ListGetLogical( BC, 'Galerkin Projector', Found )
3385+
NitscheProjector = ListGetLogical( BC,'Nitsche Projector', Found )
3386+
IF(NitscheProjector) THEN
3387+
CALL Info(Caller,'Creating an add matrix for Nitshce type of interface conditions!',Level=10)
3388+
WeakProjector = .TRUE.
3389+
ELSE
3390+
WeakProjector = ListGetLogical( BC, 'Galerkin Projector', Found )
3391+
END IF
33833392
IF (Found) THEN
33843393
StrongProjector = .NOT. WeakProjector
33853394
ELSE
@@ -3444,7 +3453,8 @@ FUNCTION LevelProjector( BMesh1, BMesh2, Repeating, AntiRepeating, &
34443453
Projector => AllocateMatrix()
34453454
Projector % FORMAT = MATRIX_LIST
34463455

3447-
IF( ListGetLogical( BC,'Nitsche Projector', Found ) ) THEN
3456+
IF( NitscheProjector ) THEN
3457+
34483458
Projector % ProjectorType = PROJECTOR_TYPE_NITSCHE
34493459
ELSE
34503460
Projector % ProjectorType = PROJECTOR_TYPE_GALERKIN
@@ -6558,9 +6568,8 @@ SUBROUTINE AddProjectorWeak1D()
65586568
! restructured into a separate routine.
65596569
IF( Projector % ProjectorType == PROJECTOR_TYPE_NITSCHE ) THEN
65606570
CALL TemporalSegmentNitscheAssembly(ElementT, NodesT, Element, Nodes, ElementM, NodesM, &
6561-
sgn0, pElemBasis, 0, &
6562-
Projector, NodeCoeff, ArcCoeff, NodeScale, NodePerm, &
6563-
InvPerm1, InvPerm2, SumArea, BC )
6571+
sgn0, pElemBasis, 0, Projector, ArcCoeff, NodeScale, NodePerm, &
6572+
InvPerm1, InvPerm2, SumArea, BC )
65646573
ELSE
65656574
CALL TemporalSegmentMortarAssembly(ElementT, NodesT, Element, Nodes, ElementM, NodesM, &
65666575
sgn0, pElemBasis, BiorthogonalBasis, CreateDual, DualMaster, DualLCoeff, 0, &

‎fem/tests/RotatingBCPoisson2D/case.sif‎

Lines changed: 2 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -5,13 +5,12 @@ Header
55
CHECK KEYWORDS Warn
66
Mesh DB "." "mortar"
77
Include Path ""
8-
Results Directory ""
8+
Results Directory "results"
99
End
1010

1111
Simulation
1212
Max Output Level = 5
1313
Coordinate System = Cartesian
14-
Coordinate Mapping(3) = 1 2 3
1514
Simulation Type = Transient
1615
Steady State Max Iterations = 1
1716
Output Intervals = 1
@@ -21,19 +20,11 @@ Simulation
2120
Timestep Sizes = 0.1
2221
Timestep Intervals = 10
2322
Output Intervals = 0
24-
! Post File = mortar.ep
25-
! Output File = case.result
23+
! Post File = case.vtu
2624

2725
Simulation Timing = Logical True
2826
End
2927

30-
Constants
31-
Gravity(4) = 0 -1 0 9.82
32-
Stefan Boltzmann = 5.67e-08
33-
Permittivity of Vacuum = 8.8542e-12
34-
Boltzmann Constant = 1.3807e-23
35-
Unit Charge = 1.602e-19
36-
End
3728

3829
Body 1
3930
Target Bodies(1) = 1
@@ -79,7 +70,6 @@ Solver 2
7970
Steady State Convergence Tolerance = 1.0e-5
8071

8172
Nonlinear System Max Iterations = 1
82-
Nonlinear System Relaxation Factor = 1.0
8373
Nonlinear System Consistent Norm = True
8474

8575
Linear System Solver = Iterative
@@ -128,15 +118,6 @@ Solver 3
128118
Parallel Operator 1 = String sum
129119
End
130120

131-
Solver 4
132-
! Exec Solver = never
133-
Exec Solver = after timestep
134-
Equation = VtuOutput
135-
Procedure = "ResultOutputSolve" "ResultOutputSolver"
136-
Output File Name = case_a
137-
Vtu Format = Logical True
138-
Single Precision = Logical True
139-
End
140121

141122
Equation 1
142123
Name = "Heat"

0 commit comments

Comments
 (0)