/[escript]/trunk/finley/src/Assemble_PDE_Single_1D.cpp
ViewVC logotype

Contents of /trunk/finley/src/Assemble_PDE_Single_1D.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6939 - (show annotations)
Mon Jan 20 03:37:18 2020 UTC (4 months, 1 week ago) by uqaeller
File size: 13796 byte(s)
Updated the copyright header.


1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2020 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Apache License, version 2.0
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014-2017 by Centre for Geoscience Computing (GeoComp)
14 * Development from 2019 by School of Earth and Environmental Sciences
15 **
16 *****************************************************************************/
17
18
19 /****************************************************************************
20
21 Assembles a single PDE into the stiffness matrix S and right hand side F
22
23 -(A_{i,j} u_,j)_i-(B_{i} u)_i+C_{j} u_,j-D u_m and -(X_,i)_i + Y
24
25 in a 1D domain. The shape functions for test and solution must be identical
26 and row_NS == row_NN.
27
28 Shape of the coefficients:
29
30 A = 1 x 1
31 B = 1
32 C = 1
33 D = scalar
34 X = 1
35 Y = scalar
36
37 *****************************************************************************/
38
39 #include "Assemble.h"
40 #include "Util.h"
41
42 #include <escript/index.h>
43
44 namespace finley {
45
46 void Assemble_PDE_Single_1D(const AssembleParameters& p,
47 const escript::Data& A, const escript::Data& B,
48 const escript::Data& C, const escript::Data& D,
49 const escript::Data& X, const escript::Data& Y)
50 {
51 const int DIM = 1;
52 bool expandedA = A.actsExpanded();
53 bool expandedB = B.actsExpanded();
54 bool expandedC = C.actsExpanded();
55 bool expandedD = D.actsExpanded();
56 bool expandedX = X.actsExpanded();
57 bool expandedY = Y.actsExpanded();
58 double *F_p = NULL;
59 if(!p.F.isEmpty()) {
60 p.F.requireWrite();
61 F_p = p.F.getSampleDataRW(0);
62 }
63 const std::vector<double>& S(p.row_jac->BasisFunctions->S);
64 const int len_EM_S = p.row_numShapesTotal*p.col_numShapesTotal;
65 const int len_EM_F = p.row_numShapesTotal;
66
67 #pragma omp parallel
68 {
69 for (index_t color = p.elements->minColor; color <= p.elements->maxColor; color++) {
70 // loop over all elements:
71 #pragma omp for
72 for (index_t e = 0; e < p.elements->numElements; e++) {
73 if (p.elements->Color[e] == color) {
74 for (int isub = 0; isub < p.numSub; isub++) {
75 const double* Vol = &p.row_jac->volume[INDEX3(0,isub,e,p.numQuadSub,p.numSub)];
76 const double* DSDX = &p.row_jac->DSDX[INDEX5(0,0,0,isub,e, p.row_numShapesTotal,DIM,p.numQuadSub,p.numSub)];
77 std::vector<double> EM_S(len_EM_S);
78 std::vector<double> EM_F(len_EM_F);
79 bool add_EM_F = false;
80 bool add_EM_S = false;
81 ///////////////
82 // process A //
83 ///////////////
84 if (!A.isEmpty()) {
85 const double *A_p = A.getSampleDataRO(e);
86 add_EM_S = true;
87 if (expandedA) {
88 const double* A_q = &(A_p[INDEX4(0,0,0,isub,DIM,DIM,p.numQuadSub)]);
89 for (int s = 0; s < p.row_numShapes; s++) {
90 for (int r = 0; r < p.col_numShapes; r++) {
91 double f = 0.;
92 for (int q = 0; q < p.numQuadSub; q++) {
93 f += Vol[q]*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*A_q[INDEX3(0,0,q,DIM,DIM)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];
94 }
95 EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=f;
96 }
97 }
98 } else { // constant A
99 for (int s = 0; s < p.row_numShapes; s++) {
100 for (int r = 0; r < p.col_numShapes; r++) {
101 double f = 0.;
102 for (int q = 0; q < p.numQuadSub; q++)
103 f += Vol[q]*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];
104 EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=f*A_p[INDEX2(0,0,DIM)];
105 }
106 }
107 }
108 }
109 ///////////////
110 // process B //
111 ///////////////
112 if (!B.isEmpty()) {
113 const double *B_p=B.getSampleDataRO(e);
114 add_EM_S=true;
115 if (expandedB) {
116 const double *B_q=&(B_p[INDEX3(0,0,isub,DIM,p.numQuadSub)]);
117 for (int s=0; s<p.row_numShapes; s++) {
118 for (int r=0; r<p.col_numShapes; r++) {
119 double f=0.;
120 for (int q=0; q<p.numQuadSub; q++) {
121 f+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*B_q[INDEX2(0,q,DIM)]*S[INDEX2(r,q,p.row_numShapes)];
122 }
123 EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=f;
124 }
125 }
126 } else { // constant B
127 for (int s=0; s<p.row_numShapes; s++) {
128 for (int r=0; r<p.col_numShapes; r++) {
129 double f=0.;
130 for (int q=0; q<p.numQuadSub; q++)
131 f+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*S[INDEX2(r,q,p.row_numShapes)];
132 EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=f*B_p[0];
133 }
134 }
135 }
136 }
137 ///////////////
138 // process C //
139 ///////////////
140 if (!C.isEmpty()) {
141 const double *C_p=C.getSampleDataRO(e);
142 add_EM_S=true;
143 if (expandedC) {
144 const double *C_q=&(C_p[INDEX3(0,0,isub,DIM,p.numQuadSub)]);
145 for (int s=0; s<p.row_numShapes; s++) {
146 for (int r=0; r<p.col_numShapes; r++) {
147 double f=0.;
148 for (int q=0; q<p.numQuadSub; q++) {
149 f+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*C_q[INDEX2(0,q,DIM)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];
150 }
151 EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=f;
152 }
153 }
154 } else { // constant C
155 for (int s=0; s<p.row_numShapes; s++) {
156 for (int r=0; r<p.col_numShapes; r++) {
157 double f=0.;
158 for (int q=0; q<p.numQuadSub; q++)
159 f+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*DSDX[INDEX3(r,0,q,p.row_numShapesTotal,DIM)];
160 EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=f*C_p[0];
161 }
162 }
163 }
164 }
165 ///////////////
166 // process D //
167 ///////////////
168 if (!D.isEmpty()) {
169 const double *D_p=D.getSampleDataRO(e);
170 add_EM_S=true;
171 if (expandedD) {
172 const double *D_q=&(D_p[INDEX2(0,isub,p.numQuadSub)]);
173 for (int s=0; s<p.row_numShapes; s++) {
174 for (int r=0; r<p.col_numShapes; r++) {
175 double f=0.;
176 for (int q=0; q<p.numQuadSub; q++) {
177 f+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*D_q[q]*S[INDEX2(r,q,p.row_numShapes)];
178 }
179 EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=f;
180 }
181 }
182 } else { // constant D
183 for (int s=0; s<p.row_numShapes; s++) {
184 for (int r=0; r<p.col_numShapes; r++) {
185 double f=0.;
186 for (int q=0; q<p.numQuadSub; q++)
187 f+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(r,q,p.row_numShapes)];
188 EM_S[INDEX4(0,0,s,r,p.numEqu,p.numComp,p.row_numShapesTotal)]+=f*D_p[0];
189 }
190 }
191 }
192 }
193 ///////////////
194 // process X //
195 ///////////////
196 if (!X.isEmpty()) {
197 const double *X_p=X.getSampleDataRO(e);
198 add_EM_F=true;
199 if (expandedX) {
200 const double *X_q=&(X_p[INDEX3(0,0,isub,DIM,p.numQuadSub)]);
201 for (int s=0; s<p.row_numShapes; s++) {
202 double f=0.;
203 for (int q=0; q<p.numQuadSub; q++)
204 f+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)]*X_q[INDEX2(0,q,DIM)];
205 EM_F[INDEX2(0,s,p.numEqu)]+=f;
206 }
207 } else { // constant X
208 for (int s=0; s<p.row_numShapes; s++) {
209 double f=0.;
210 for (int q=0; q<p.numQuadSub; q++)
211 f+=Vol[q]*DSDX[INDEX3(s,0,q,p.row_numShapesTotal,DIM)];
212 EM_F[INDEX2(0,s,p.numEqu)]+=f*X_p[0];
213 }
214 }
215 }
216 ///////////////
217 // process Y //
218 ///////////////
219 if (!Y.isEmpty()) {
220 const double *Y_p=Y.getSampleDataRO(e);
221 add_EM_F=true;
222 if (expandedY) {
223 const double *Y_q=&(Y_p[INDEX2(0,isub, p.numQuadSub)]);
224 for (int s = 0; s < p.row_numShapes; s++) {
225 double f = 0.;
226 for (int q = 0; q < p.numQuadSub; q++)
227 f += Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*Y_q[q];
228 EM_F[INDEX2(0,s,p.numEqu)]+=f;
229 }
230 } else { // constant Y
231 for (int s = 0; s < p.row_numShapes; s++) {
232 double f = 0.;
233 for (int q = 0; q < p.numQuadSub; q++)
234 f+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];
235 EM_F[INDEX2(0,s,p.numEqu)]+=f*Y_p[0];
236 }
237 }
238 }
239 // add the element matrices onto the matrix and
240 // right hand side
241 IndexVector row_index(p.row_numShapesTotal);
242 for (int q = 0; q < p.row_numShapesTotal; q++)
243 row_index[q] = p.row_DOF[p.elements->Nodes[INDEX2(p.row_node[INDEX2(q, isub, p.row_numShapesTotal)], e, p.NN)]];
244
245 if (add_EM_F)
246 util::addScatter(p.row_numShapesTotal,
247 &row_index[0], p.numEqu, &EM_F[0], F_p,
248 p.row_DOF_UpperBound);
249 if (add_EM_S)
250 Assemble_addToSystemMatrix(p.S,
251 p.row_numShapesTotal, &row_index[0],
252 p.numEqu, p.col_numShapesTotal,
253 &row_index[0], p.numComp, &EM_S[0]);
254 } // end of isub
255 } // end color check
256 } // end element loop
257 } // end color loop
258 } // end parallel region
259 }
260
261 } // namespace finley
262

Properties

Name Value
svn:mergeinfo /branches/4.0fordebian/finley/src/Assemble_PDE_Single_1D.cpp:5567-5588 /branches/lapack2681/finley/src/Assemble_PDE_Single2_1D.cpp:2682-2741 /branches/pasowrap/finley/src/Assemble_PDE_Single2_1D.cpp:3661-3674 /branches/py3_attempt2/finley/src/Assemble_PDE_Single2_1D.cpp:3871-3891 /branches/restext/finley/src/Assemble_PDE_Single2_1D.cpp:2610-2624 /branches/ripleygmg_from_3668/finley/src/Assemble_PDE_Single2_1D.cpp:3669-3791 /branches/stage3.0/finley/src/Assemble_PDE_Single2_1D.cpp:2569-2590 /branches/symbolic_from_3470/finley/src/Assemble_PDE_Single2_1D.cpp:3471-3974 /branches/trilinos_from_5897/finley/src/Assemble_PDE_Single_1D.cpp:5898-6118 /release/3.0/finley/src/Assemble_PDE_Single2_1D.cpp:2591-2601 /release/4.0/finley/src/Assemble_PDE_Single_1D.cpp:5380-5406 /trunk/finley/src/Assemble_PDE_Single2_1D.cpp:4257-4344

  ViewVC Help
Powered by ViewVC 1.1.26