CoCalc -- PoissonAssembly.f90

GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/fem/tests/BlockPoisson2/PoissonAssembly.f90
³²⁰⁶ views
1

2
!------------------------------------------------------------------------------
3
SUBROUTINE BulkAssembly( Model,Solver,dt,Transient, &
4
    Mass, Damp, Stiff, Force, InElement, nrow, ncol )
5
!------------------------------------------------------------------------------
6

7
   USE DefUtils
8

9
 
10
  IMPLICIT NONE
11
!------------------------------------------------------------------------------
12
  TYPE(Model_t)  :: Model
13
  TYPE(Solver_t), TARGET :: Solver
14
  LOGICAL ::  Transient
15
  REAL(KIND=dp) :: dt
16
  REAL(KIND=dp) :: Stiff(:,:), Damp(:,:), Mass(:,:), Force(:)
17
  TYPE(Element_t), TARGET :: InElement
18
  INTEGER :: Nrow, Ncol 
19
!------------------------------------------------------------------------------
20
!    Local variables
21
!------------------------------------------------------------------------------
22

23
  TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff
24
  TYPE(Nodes_t) :: Nodes
25
  REAL(KIND=dp), ALLOCATABLE :: Source(:), Cond(:), React(:)
26
  REAL(KIND=dp), ALLOCATABLE :: Basis(:), dBasisdx(:,:)
27
  REAL(KIND=dp) :: SourceAtIP, CondAtIP, ReactAtIp, Weight, DetJ
28
  INTEGER :: i,j,t,p,q,n
29
  LOGICAL :: Visited = .FALSE., Found
30
  TYPE(Element_t), POINTER :: Element
31
  TYPE(ValueList_t), POINTER :: Material
32

33
  
34
  SAVE Visited, Source, Cond, React, Nodes, Basis, dBasisdx
35

36

37
  IF( .NOT. Visited ) THEN
38
     N = Solver % Mesh % MaxElementNodes    
39
     ALLOCATE( Source(n), Cond(n), React(n), Basis(n), dBasisdx(n,3), &
40
         Nodes % x(n), Nodes % y(n), Nodes % z(n) )     
41
     CALL info('PoissonBulkAssembly','1st time')
42
     Visited = .TRUE.
43
  END IF
44

45
  n = Nrow 
46
  CALL GetElementNodes( Nodes ) 
47
  Material => GetMaterial()  
48
  React(1:n) = GetReal( Material,'Reaction',Found )
49
  Source(1:n) = GetReal( Material,'Source' )
50
  Cond(1:n) = GetReal( Material,'Conductivity')
51

52
  Element => GetCurrentElement()
53
  IntegStuff = GaussPoints( Element )
54

55
  DO t=1,IntegStuff % n
56
    
57
    Found = ElementInfo( Element, Nodes, IntegStuff % u(t), &
58
        IntegStuff % v(t), IntegStuff % w(t), detJ, Basis, dBasisdx )
59

60
    Weight = IntegStuff % s(t) * detJ
61

62
    CondAtIP = SUM( Basis(1:Nrow) * Cond(1:Nrow) )
63
    SourceAtIP = SUM( Basis(1:Nrow) * Source(1:Nrow) )
64
    ReactAtIP = SUM( Basis(1:Nrow) * React(1:Nrow) )
65
    
66
    DO p=1,n
67
      DO q=1,n
68
        STIFF(p,q) = STIFF(p,q) + Weight * CondAtIP * &
69
            SUM(dBasisdx(p,1:3)*dBasisdx(q,1:3))
70
        STIFF(p,q) = STIFF(p,q) + Weight * ReactAtIP * &
71
            Basis(p)*Basis(q)
72
      END DO
73
      FORCE(p) = FORCE(p) + Weight * SourceAtIP * Basis(p)
74
    END DO
75
  END DO
76
  
77
!------------------------------------------------------------------------------
78
END SUBROUTINE BulkAssembly
79
!------------------------------------------------------------------------------
80

81

82

83
!------------------------------------------------------------------------------
84
SUBROUTINE BoundaryAssembly( Model,Solver,dt,Transient, &
85
            Mass, Damp, Stiff, Force, InElement, nrow, ncol )
86
!------------------------------------------------------------------------------
87

88
  USE DefUtils
89
 
90
  IMPLICIT NONE
91
!------------------------------------------------------------------------------
92
  TYPE(Model_t)  :: Model
93
  TYPE(Solver_t), TARGET :: Solver
94
  LOGICAL ::  Transient
95
  REAL(KIND=dp) :: dt
96
  REAL(KIND=dp) :: Stiff(:,:), Damp(:,:), Mass(:,:), Force(:)
97
  TYPE(Element_t), TARGET :: InElement
98
  INTEGER :: Nrow, Ncol 
99
!------------------------------------------------------------------------------
100
!    Local variables
101
!------------------------------------------------------------------------------
102

103
  TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff
104
  TYPE(Nodes_t) :: Nodes
105
  REAL(KIND=dp), ALLOCATABLE :: Flux(:), Text(:), Coeff(:)
106
  REAL(KIND=dp), ALLOCATABLE :: Basis(:), dBasisdx(:,:)
107
  REAL(KIND=dp) :: FluxAtIP, TextAtIP, CoeffAtIP, Weight, DetJ
108
  INTEGER :: i,j,t,p,q,n
109
  LOGICAL :: Visited = .FALSE., Found, Found2
110
  TYPE(ValueList_t), POINTER :: BC
111
  TYPE(Element_t), POINTER :: Element
112

113
  
114
  SAVE Visited, Flux, Text, Coeff, Nodes, Basis, dBasisdx
115

116

117
  IF( .NOT. Visited ) THEN
118
    N = Solver % Mesh % MaxElementNodes    
119
    ALLOCATE( Flux(n), Text(n), Coeff(n), Basis(n), dBasisdx(n,3), &
120
        Nodes % x(n), Nodes % y(n), Nodes % z(n) )     
121
    CALL info('PoissonBoundaryAssembly','1st time')
122
    Visited = .TRUE.
123
  END IF
124
  
125
  Element => GetCurrentElement()
126
  n = Nrow 
127
  CALL GetElementNodes( Nodes ) 
128
  
129
  BC => GetBC( )
130
  
131
  Flux(1:n) = GetReal( BC,'Flux',Found )
132
  Coeff(1:n) = GetReal( BC,'Coefficient',Found2 )
133
  IF(.NOT. (Found .OR. Found2)) RETURN
134
  Text(1:n) = GetReal( BC,'External',Found )
135

136

137
  IntegStuff = GaussPoints( Element )
138
  
139
  DO t=1,IntegStuff % n
140
    
141
    Found = ElementInfo( Element, Nodes, IntegStuff % u(t), &
142
        IntegStuff % v(t), IntegStuff % w(t), detJ, Basis, dBasisdx )
143
    
144
    Weight = IntegStuff % s(t) * detJ
145
    
146
    FluxAtIP = SUM( Basis(1:Nrow) * Flux(1:Nrow) )
147
    CoeffAtIP = SUM( Basis(1:Nrow) * Coeff(1:Nrow) )
148
    TextAtIP = SUM( Basis(1:Nrow) * Text(1:Nrow) )
149
    
150
    DO p=1,n
151
      DO q=1,n
152
        STIFF(p,q) = STIFF(p,q) + Weight * CoeffAtIP * &
153
            Basis(p)*Basis(q)
154
      END DO
155
      FORCE(p) = FORCE(p) + Weight * Basis(p) * &
156
          ( FluxAtIP + CoeffAtIP * TextAtIP )
157
    END DO
158
  END DO
159

160
!------------------------------------------------------------------------------
161
END SUBROUTINE BoundaryAssembly
162
!------------------------------------------------------------------------------
163

164

165
Product

Resources

Company
