1SUBROUTINE EdgeElementSolver( Model,Solver,dt,TransientSimulation )
2!------------------------------------------------------------------------------
3!******************************************************************************
4!
5!  Solve the best approximation of the vector field U = (1+z-y,1-z+x,1-x+y)
6!  with respect to the L2 norm using H(curl)-conforming basis functions.
7!  Additionally, compute the relative error of the solution using
8!  the energy norm corresponding to the operator I - curl curl.
9!  For the basis functions obtained via the function EdgeElementInfo the exact
10!  solution should be in the FE solution space. This solver thus offers
11!  a consistency check for creating discretizations based on the basis functions
12!  provided by the function EdgeElementInfo. This test can be performed basically
13!  on any 3-D mesh if the element definition of the sif file is adjusted accordingly.
14!
15!  ARGUMENTS:
16!
17!  TYPE(Model_t) :: Model,
18!     INPUT: All model information (mesh, materials, BCs, etc...)
19!
20!  TYPE(Solver_t) :: Solver
21!     INPUT: Linear & nonlinear equation solver options
22!
23!  REAL(KIND=dp) :: dt
24!     INPUT: Timestep size for time dependent simulations
25!
26!  LOGICAL :: TransientSimulation
27!     INPUT: Steady state or transient simulation
28!
29!******************************************************************************
30  USE DefUtils
31
32  IMPLICIT NONE
33!------------------------------------------------------------------------------
34  TYPE(Solver_t) :: Solver
35  TYPE(Model_t) :: Model
36
37  REAL(KIND=dp) :: dt
38  LOGICAL :: TransientSimulation
39!------------------------------------------------------------------------------
40! Local variables
41!------------------------------------------------------------------------------
42  LOGICAL :: AllocationsDone = .FALSE., Found
43  TYPE(Element_t), POINTER :: Element
44  TYPE(Nodes_t) :: Nodes
45  TYPE(ValueList_t), POINTER :: BodyForce, Material, BC
46  TYPE(Variable_t), POINTER :: Var
47
48  REAL(KIND=dp) :: Norm, u, v, w, Err, EK, SolNorm
49
50  INTEGER :: n, ne, nf, nb, np, nd, t, istat, i, j, k, l, active, dim
51
52  REAL(KIND=dp), ALLOCATABLE :: LOAD(:,:), Acoef(:)
53  REAL(KIND=dp), ALLOCATABLE :: STIFF(:,:), FORCE(:)
54
55  TYPE(Mesh_t), POINTER :: Mesh
56  TYPE(Matrix_t), POINTER :: A
57
58  LOGICAL :: stat, PiolaVersion, ErrorEstimation, UseTabulatedBasis
59  LOGICAL :: UseEnergyNorm
60
61  INTEGER, ALLOCATABLE :: Indices(:)
62
63  SAVE STIFF, LOAD, FORCE, Acoef, AllocationsDone, Nodes, Indices
64!------------------------------------------------------------------------------
65  PiolaVersion = .TRUE.
66  dim = CoordinateSystemDimension()
67
68  !Allocate some permanent storage, this is done first time only:
69  !--------------------------------------------------------------
70  Mesh => GetMesh()
71
72  IF ( .NOT. AllocationsDone ) THEN
73     N = Mesh % MaxElementDOFs  ! just big enough
74     ALLOCATE( FORCE(N), LOAD(6,N), STIFF(N,N), &
75          Acoef(N), Indices(N), STAT=istat )
76     IF ( istat /= 0 ) THEN
77        CALL Fatal( 'EdgeElementSolver', 'Memory allocation error.' )
78     END IF
79     AllocationsDone = .TRUE.
80  END IF
81
82  Solver % Matrix % COMPLEX = .FALSE.
83  A => GetMatrix()
84
85  ErrorEstimation = GetLogical( GetSolverParams(), 'Error Computation', Found)
86  UseEnergyNorm = GetLogical( GetSolverParams(), 'Use Energy Norm', Found)
87  UseTabulatedBasis = GetLogical( GetSolverParams(), 'Tabulate Basis', Found)
88
89  !-----------------------
90  ! System assembly:
91  !----------------------
92  active = GetNOFActive()
93  CALL DefaultInitialize()
94
95  DO t=1,active
96     Element => GetActiveElement(t)
97     n  = GetElementNOFNodes() ! The number of nodes corresponding to the background mesh
98     nd = GetElementNOFDOFs()
99     nb = GetElementNOFBDOFs()
100
101     LOAD = 0.0d0
102     IF (.FALSE.) THEN
103        BodyForce => GetBodyForce()
104        IF ( ASSOCIATED(BodyForce) ) THEN
105           Load(1,1:n) = GetReal( BodyForce, 'Body Source 1', Found )
106           Load(2,1:n) = GetReal( BodyForce, 'Body Source 2', Found )
107           Load(3,1:n) = GetReal( BodyForce, 'Body Source 3', Found )
108        END IF
109     END IF
110
111     Acoef(1:n) = 1.0d0
112     IF (.FALSE.) THEN
113        Material => GetMaterial( Element )
114        IF ( ASSOCIATED(Material) ) THEN
115           Acoef(1:n) = GetReal( Material, 'Conductivity', Found )
116           IF (.NOT. Found) CALL Fatal( 'NedelecSolve', 'Conductivity must be specified' )
117        END IF
118     END IF
119
120     ! Perform an additional check that DOF counts are right:
121     !-------------------------------------------------------------------
122     IF (GetElementFamily() == 5 .AND. nd /= 6) THEN
123        WRITE(Message,'(I2,A)') nd, 'DOFs Found'
124        CALL Fatal('EdgeElementSolver','Indices for a tetrahedron erratic')
125     END IF
126     IF (GetElementFamily() == 6 .AND. nd /= 10) THEN
127        WRITE(Message,'(I2,A)') nd, 'DOFs Found'
128        CALL Fatal('EdgeElementSolver','Indices for a pyramid erratic')
129     END IF
130     IF (GetElementFamily() == 7 .AND. nd /= 15) THEN
131        WRITE(Message,'(I2,A)') nd, 'DOFs Found'
132        CALL Fatal('EdgeElementSolver','Indices for a prism erratic')
133     END IF
134     IF (GetElementFamily() == 8 .AND. nd /= 27) THEN
135       WRITE(Message,'(I2,A)') nd, 'DOFs Found'
136       CALL Fatal('EdgeElementSolver','Indices for a brick erratic')
137     END IF
138
139
140     !Get element local matrix and rhs vector:
141     !----------------------------------------
142     CALL LocalMatrix( STIFF, FORCE, LOAD, Acoef, Element, n, nd+nb, dim, UseTabulatedBasis)
143
144     !Update global matrix and rhs vector from local matrix & vector:
145     !---------------------------------------------------------------
146     CALL DefaultUpdateEquations( STIFF, FORCE )
147
148  END DO
149
150  CALL DefaultFinishBulkAssembly()
151  CALL DefaultFinishAssembly()
152
153  CALL DefaultDirichletBCs()
154
155  Norm = DefaultSolve()
156
157  !-------------------------------------------------------------------
158  ! Compute the error norm for the model problem considered. Note that
159  ! this part has not been modified to support the option
160  ! Tabulate Basis = Logical True.
161  !--------------------------------------------------------------------
162  IF (ErrorEstimation) THEN
163     Err = 0.0d0
164     SolNorm = 0.0d0
165     DO t=1,Solver % NumberOfActiveElements
166        Element => GetActiveElement(t)
167        n  = GetElementNOFNodes()
168        nd = GetElementDOFs( Indices )
169
170        Load(1,1:nd) = Solver % Variable % Values( Solver % Variable % &
171             Perm(Indices(1:nd)) )
172
173        CALL MyComputeError(Load, Element, n, nd, dim, Err, SolNorm, UseEnergyNorm)
174     END DO
175
176     WRITE (*, '(A,E16.8)') 'Error Norm = ', SQRT(ParallelReduction(Err))/SQRT(ParallelReduction(SolNorm))
177  END IF
178
179CONTAINS
180
181
182
183
184!---------------------------------------------------------------------------------------------
185  SUBROUTINE LocalMatrix(  STIFF, FORCE, LOAD, Acoef, Element, n, nd, dim, UseTabulatedBasis)
186!---------------------------------------------------------------------------------------------
187    REAL(KIND=dp) :: STIFF(:,:), FORCE(:), LOAD(:,:), Acoef(:)
188    TYPE(Element_t), POINTER :: Element
189    INTEGER :: n, nd, dim
190    LOGICAL :: UseTabulatedBasis
191    !------------------------------------------------------------------------------
192    REAL(KIND=dp) :: nu
193    REAL(KIND=dp) :: EBasis(nd,3), CurlEBasis(nd,3), F(3,3), G(3,3)
194    REAL(KIND=dp) :: dBasisdx(n,3)
195    REAL(KIND=dp) :: Basis(n), DetJ, xq, yq, zq, uq, vq, wq, sq
196    LOGICAL :: Stat, Found
197    INTEGER :: t, i, j, p, q, np
198
199    TYPE(GaussIntegrationPoints_t) :: IP
200    TYPE(Nodes_t), SAVE :: Nodes
201    ! ----------------------------------------------------------------------------
202    INTEGER :: PermVec(nd)
203    REAL(KIND=dp) :: SignVec(nd)
204    REAL(KIND=dp) :: ReadyBasis(nd,3), ReadyCurlBasis(nd,3)
205
206    LOGICAL, SAVE :: BricksVisited = .FALSE., PrismsVisited = .FALSE.
207    LOGICAL, SAVE :: PyramidsVisited = .FALSE., TetraVisited = .FALSE.
208
209    REAL(KIND=dp), TARGET :: TetraBasis(6,12), TetraCurlBasis(6,12)
210    REAL(KIND=dp), TARGET :: PyramidBasis(10,36), PyramidCurlBasis(10,36)
211    REAL(KIND=dp), TARGET :: PrismBasis(15,54), PrismCurlBasis(15,54)
212    REAL(KIND=dp), TARGET :: BrickBasis(27,81), BrickCurlBasis(27,81)
213    REAL(KIND=dp), POINTER :: BasisTable(:,:), CurlBasisTable(:,:)
214
215    SAVE PrismBasis, PrismCurlBasis, BrickBasis, BrickCurlBasis
216    SAVE TetraBasis, TetraCurlBasis, PyramidBasis, PyramidCurlBasis
217    !------------------------------------------------------------------------------
218    CALL GetElementNodes( Nodes )
219
220    STIFF = 0.0d0
221    FORCE = 0.0d0
222
223    !-------------------------------------
224    ! Numerical integration over element:
225    !-------------------------------------
226    IP = GaussPoints(Element, EdgeBasis=.TRUE., PReferenceElement=PiolaVersion)
227
228    IF (UseTabulatedBasis .AND. PiolaVersion) THEN
229       !----------------------------------------------------------------------------
230       ! Obtain the data for permuting basis function positions and applying
231       ! sign reversions. This data is the same for all integration points.
232       ! If elements of this type has not yet been visited, tabulate basis
233       ! function values to avoid the recomputation.
234       !----------------------------------------------------------------------------
235       IF ( (GetElementFamily() == 5 .AND. TetraVisited) .OR. &
236            (GetElementFamily() == 6 .AND. PyramidsVisited) .OR. &
237            (GetElementFamily() == 7 .AND. PrismsVisited) .OR. &
238            (GetElementFamily() == 8 .AND. BricksVisited) ) THEN
239          CALL ReorderingAndSignReversionsData(Element,Nodes,PermVec,SignVec)
240
241       ELSE
242          !---------------------------------------------------------------------------
243          ! The first visit for this element type, tabulate the basis function values.
244          !---------------------------------------------------------------------------
245          CALL ReorderingAndSignReversionsData(Element,Nodes,PermVec,SignVec)
246          DO t=1,IP % n
247             stat = EdgeElementInfo( Element, Nodes, IP % U(t), IP % V(t), &
248                  IP % W(t), DetF=detJ, Basis=Basis, EdgeBasis=EBasis, RotBasis=CurlEBasis, &
249                  ApplyPiolaTransform = .FALSE.)
250             !------------------------------------------------------------------------
251             ! Revert order and sign changes to the reference element default
252             ! and tabulate the values for later usage
253             !------------------------------------------------------------------------
254             DO i=1,nd
255                EBasis(i,1:3) = SignVec(i) * EBasis(i,1:3)
256                CurlEBasis(i,1:3) = SignVec(i) * CurlEBasis(i,1:3)
257             END DO
258
259             SELECT CASE( GetElementFamily() )
260             CASE(5)
261                DO i=1,nd
262                   j = PermVec(i)
263                   TetraBasis(i,(t-1)*3+1:(t-1)*3+3) = EBasis(j,1:3)
264                   TetraCurlBasis(i,(t-1)*3+1:(t-1)*3+3) = CurlEBasis(j,1:3)
265                END DO
266             CASE(6)
267                DO i=1,nd
268                   j = PermVec(i)
269                   PyramidBasis(i,(t-1)*3+1:(t-1)*3+3) = EBasis(j,1:3)
270                   PyramidCurlBasis(i,(t-1)*3+1:(t-1)*3+3) = CurlEBasis(j,1:3)
271                END DO
272             CASE(7)
273                DO i=1,nd
274                   j = PermVec(i)
275                   PrismBasis(i,(t-1)*3+1:(t-1)*3+3) = EBasis(j,1:3)
276                   PrismCurlBasis(i,(t-1)*3+1:(t-1)*3+3) = CurlEBasis(j,1:3)
277                END DO
278             CASE(8)
279                DO i=1,nd
280                   j = PermVec(i)
281                   BrickBasis(i,(t-1)*3+1:(t-1)*3+3) = EBasis(j,1:3)
282                   BrickCurlBasis(i,(t-1)*3+1:(t-1)*3+3) = CurlEBasis(j,1:3)
283                END DO
284             END SELECT
285
286          END DO
287
288          SELECT CASE( GetElementFamily() )
289          CASE(5)
290             TetraVisited = .true.
291             CALL Info( 'EdgeElementSolver', 'TETRAHEDRAL BASIS FUNCTIONS WERE TABULATED', Level=10)
292             WRITE(Message,'(A,I2)') 'Integration points ', IP % n
293             CALL Info('EdgeElementSolver', Message, Level=10)
294          CASE(6)
295             PyramidsVisited = .TRUE.
296             CALL Info( 'EdgeElementSolve', 'PYRAMIDICAL BASIS FUNCTIONS WERE TABULATED', Level=10)
297             WRITE(Message,'(A,I2)') 'Integration points ', IP % n
298             CALL Info('EdgeElementSolver', Message, Level=10)
299          CASE(7)
300             PrismsVisited = .TRUE.
301             CALL Info( 'EdgeElementSolve', 'PRISM BASIS FUNCTIONS WERE TABULATED', Level=10)
302             WRITE(Message,'(A,I2)') 'Integration points ', IP % n
303             CALL Info('EdgeElementSolver', Message, Level=10)
304          CASE(8)
305             BricksVisited = .TRUE.
306             CALL Info( 'EdgeElementSolve', 'BRICK BASIS FUNCTIONS WERE TABULATED', Level=10)
307             WRITE(Message,'(A,I2)') 'Integration points ', IP % n
308             CALL Info('EdgeElementSolver', Message, Level=10)
309          END SELECT
310       END IF
311    END IF
312
313    np = 0  ! Set np = n, if nodal dofs are employed; otherwise set np = 0
314
315    DO t=1,IP % n
316
317       IF (PiolaVersion) THEN
318          IF (UseTabulatedBasis) THEN
319             SELECT CASE(Element % TYPE % ElementCode / 100)
320             CASE(5)
321                BasisTable => TetraBasis(1:6,(t-1)*3+1:(t-1)*3+3)
322                CurlBasisTable => TetraCurlBasis(1:6,(t-1)*3+1:(t-1)*3+3)
323             CASE(6)
324                BasisTable => PyramidBasis(1:10,(t-1)*3+1:(t-1)*3+3)
325                CurlBasisTable => PyramidCurlBasis(1:10,(t-1)*3+1:(t-1)*3+3)
326             CASE(7)
327                BasisTable => PrismBasis(1:15,(t-1)*3+1:(t-1)*3+3)
328                CurlBasisTable => PrismCurlBasis(1:15,(t-1)*3+1:(t-1)*3+3)
329             CASE(8)
330                BasisTable => BrickBasis(1:27,(t-1)*3+1:(t-1)*3+3)
331                CurlBasisTable => BrickCurlBasis(1:27,(t-1)*3+1:(t-1)*3+3)
332             CASE DEFAULT
333                CALL Fatal( 'EdgeElementSolver', 'THE BASIS FUNCTIONS FOR THIS ELEMENT TYPE NONTABULATED' )
334             END SELECT
335             !----------------------------------------------------------------
336             ! Permute, apply sign reversions and apply the Piola transform via
337             ! calling EdgeElementInfo:
338             !----------------------------------------------------------------
339             DO i=1,nd
340                ReadyBasis(PermVec(i),1:3) = BasisTable(i,1:3)
341                ReadyCurlBasis(PermVec(i),1:3) = CurlBasisTable(i,1:3)
342             END DO
343             DO i=1,nd
344                ReadyBasis(i,1:3) = SignVec(i) * ReadyBasis(i,1:3)
345                ReadyCurlBasis(i,1:3) = SignVec(i) * ReadyCurlBasis(i,1:3)
346             END DO
347             stat = EdgeElementInfo( Element, Nodes, IP % U(t), IP % V(t), &
348                  IP % W(t), DetF=detJ, Basis=Basis, EdgeBasis=EBasis, &
349                  RotBasis=CurlEBasis, ApplyPiolaTransform = .TRUE., &
350                  ReadyEdgeBasis=ReadyBasis, ReadyRotBasis = ReadyCurlBasis)
351
352          ELSE
353             stat = EdgeElementInfo( Element, Nodes, IP % U(t), IP % V(t), &
354                  IP % W(t), detF=detJ, Basis=Basis, EdgeBasis=EBasis, &
355                  RotBasis=CurlEBasis, ApplyPiolaTransform = .TRUE.)
356          END IF
357
358       ELSE
359          stat = ElementInfo( Element, Nodes, IP % U(t), IP % V(t), &
360               IP % W(t), detJ, Basis, dBasisdx )
361          CALL GetEdgeBasis(Element, EBasis, CurlEBasis, Basis, dBasisdx)
362       END IF
363
364       xq = SUM( Nodes % x(1:n) * Basis(1:n) )
365       yq = SUM( Nodes % y(1:n) * Basis(1:n) )
366       zq = SUM( Nodes % z(1:n) * Basis(1:n) )
367
368       !nu = SUM( Basis(1:n) * Acoef(1:n) )
369
370       !----------------------------------------------------------------
371       ! The following branch could be used to produce the
372       ! Galerkin projection of a solution component for visualization.
373       !------------------------------------------------------------------
374       IF (np > 0) THEN
375          DO p = 1,n
376             DO q = 1,n
377                STIFF(p,q) = STIFF(p,q) + Basis(p) * Basis(q) * detJ * IP % s(t)
378             END DO
379
380             DO q = 1,nd-np
381                j = np + q
382                ! The following is for plotting the x-component of the solution
383                STIFF(p,j) = STIFF(p,j) - EBasis(q,1) * Basis(p) * detJ * IP % s(t)
384             END DO
385          END DO
386       END IF
387
388       !--------------------------------------------------------------
389       ! The equations for H(curl)-conforming part...
390       !---------------------------------------------------------------
391       DO p = 1,nd-np
392          !----------------------------
393          ! The inner product (u,v)_K
394          !----------------------------
395          i = np + p
396          DO q = 1,nd-np
397             j = np + q
398             STIFF(i,j) = STIFF(i,j) + 1.0d0 * &
399                  SUM( EBasis(q,1:dim) * EBasis(p,1:dim) ) * detJ * IP % s(t)
400          END DO
401
402          !----------------------------------------
403          ! RHS corresponding to the exact solution
404          !----------------------------------------
405          FORCE(i) = FORCE(i) +  1.0d0 * EBasis(p,1) * detJ * IP % s(t) + &
406               1.0d0 * EBasis(p,2) * detJ * IP % s(t) + &
407               1.0d0 * EBasis(p,3) * detJ * IP % s(t) + &
408               (zq * EBasis(p,1) - xq *  EBasis(p,3)) * detJ * IP % s(t) + &
409               (yq * EBasis(p,3) - zq *  EBasis(p,2)) * detJ * IP % s(t) + &
410               (-yq * EBasis(p,1) + xq *  EBasis(p,2)) * detJ * IP % s(t)
411       END DO
412    END DO
413!------------------------------------------------------------------------------
414  END SUBROUTINE LocalMatrix
415!------------------------------------------------------------------------------
416
417
418
419!----------------------------------------------------------------------------------
420  SUBROUTINE MyComputeError(LOAD, Element, n, nd, dim, EK, SolNorm, UseEnergyNorm)
421!----------------------------------------------------------------------------------
422    REAL(KIND=dp) :: Load(:,:), EK, SolNorm
423    TYPE(Element_t), POINTER :: Element
424    INTEGER :: n, nd, dim
425    LOGICAL :: UseEnergyNorm
426!--------------------------------------------------------------------------------
427    REAL(KIND=dp) :: EBasis(nd,3), CurlEBasis(nd,3)
428    REAL(KIND=dp) :: Basis(n), DetJ, xq, yq, zq, uq, vq, wq, sq, &
429         u(3), rotu(3), sol(3), rotsol(3), e(3), rote(3), F(3,3), G(3,3)
430    REAL(KIND=dp) :: dBasisdx(n,3)
431    LOGICAL :: Stat, Found
432    INTEGER :: t, i, j, p, q, np
433    TYPE(GaussIntegrationPoints_t) :: IP
434
435    TYPE(Nodes_t), SAVE :: Nodes
436    !------------------------------------------------------------------------------
437    CALL GetElementNodes( Nodes )
438
439    !-------------------------------------
440    ! Numerical integration over element:
441    !-------------------------------------
442    IP = GaussPoints(Element, EdgeBasis=.TRUE., PReferenceElement=PiolaVersion)
443
444    np = 0  ! Set np = n, if nodal dofs are employed; otherwise set np = 0
445
446    DO t=1,IP % n
447
448       IF (PiolaVersion) THEN
449          stat = EdgeElementInfo( Element, Nodes, IP % U(t), IP % V(t), &
450               IP % W(t), F, G, detJ, Basis, EBasis, CurlEBasis, ApplyPiolaTransform = .TRUE.)
451       ELSE
452          stat = ElementInfo( Element, Nodes, IP % U(t), IP % V(t), &
453               IP % W(t), detJ, Basis, dBasisdx )
454          CALL GetEdgeBasis(Element, EBasis, CurlEBasis, Basis, dBasisdx)
455       END IF
456
457       xq = SUM( Nodes % x(1:n) * Basis(1:n) )
458       yq = SUM( Nodes % y(1:n) * Basis(1:n) )
459       zq = SUM( Nodes % z(1:n) * Basis(1:n) )
460
461       u(1) = SUM( Load(1,np+1:nd) * EBasis(1:nd-np,1) )
462       u(2) = SUM( Load(1,np+1:nd) * EBasis(1:nd-np,2) )
463       u(3) = SUM( Load(1,np+1:nd) * EBasis(1:nd-np,3) )
464
465       rotu(1) = SUM( Load(1,np+1:nd) * CurlEBasis(1:nd-np,1) )
466       rotu(2) = SUM( Load(1,np+1:nd) * CurlEBasis(1:nd-np,2) )
467       rotu(3) = SUM( Load(1,np+1:nd) * CurlEBasis(1:nd-np,3) )
468
469       ! Compute the square of the energy norm of the solution and error:
470       sol = 0.0d0
471       sol(1) = 1.0d0 + zq - yq
472       sol(2) = 1.0d0 - zq + xq
473       sol(3) = 1.0d0 - xq + yq
474       rotsol(1:3) = 2.0d0
475       e(:) = sol(:) - u(:)
476       rote(:) = rotsol(:) - rotu(:)
477
478       IF (UseEnergyNorm) THEN
479          ! Energy norm
480          SolNorm = SolNorm + (SUM( Sol(1:3) * Sol(1:3) ) + SUM( rotsol(1:3) * rotsol(1:3) )) * detJ
481          EK = EK + (SUM( e(1:3) * e(1:3) ) + SUM( rote(1:3) * rote(1:3) )) * detJ
482       ELSE
483          ! L2 norm
484          SolNorm = SolNorm + SUM( Sol(1:3) * Sol(1:3) )* detJ
485          EK = EK + SUM( e(1:3) * e(1:3) )* detJ
486       END IF
487    END DO
488!------------------------------------------------------------------------------
489  END SUBROUTINE MyComputeError
490!-----------------------------------------------------------------------------
491
492!------------------------------------------------------------------------------
493END SUBROUTINE EdgeElementSolver
494!------------------------------------------------------------------------------
495