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

Contents of /trunk/finley/src/Assemble_PDE.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: 12165 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 Assembles the system of numEqu PDEs into the stiffness matrix S and right
21 hand side F:
22
23 -div(A*grad u)-div(B*u)+C*grad u + D*u = -div X + Y
24
25 -(A_{k,i,m,j} u_m,j)_i-(B_{k,i,m} u_m)_i+C_{k,m,j} u_m,j-D_{k,m} u_m = -(X_{k,i})_i + Y_k
26
27 u has numComp components.
28 Shape of the coefficients:
29
30 A = numEqu x numDim x numComp x numDim
31 B = numDim x numEqu x numComp
32 C = numEqu x numDim x numComp
33 D = numEqu x numComp
34 X = numEqu x numDim
35 Y = numEqu
36
37 The coefficients A,B,C,D,X and Y have to be defined on the integration
38 points or not present (i.e. empty).
39
40 S and F have to be initialized before the routine is called. S or F can
41 be NULL. In this case the left or the right hand side of the PDE is not
42 processed.
43
44 The routine does not consider any boundary conditions.
45
46 *****************************************************************************/
47
48 #include "Assemble.h"
49 #include "Util.h"
50
51 #include <sstream>
52
53 namespace finley {
54
55 using escript::DataTypes::real_t;
56 using escript::DataTypes::cplx_t;
57
58 inline void setNumSamplesError(const char* c, int n0, int n1)
59 {
60 std::stringstream ss;
61 ss << "Assemble_PDE: number of sample points of coefficient " << c
62 << " don't match (" << n0 << "," << n1 << ").";
63 const std::string errorMsg(ss.str());
64 throw escript::ValueError(errorMsg);
65 }
66
67 inline void setShapeError(const char* c, int num, const int* dims)
68 {
69 std::stringstream ss;
70 ss << "Assemble_PDE: shape of coefficient " << c
71 << " does not match (" << dims[0] << ",";
72 if (num > 1) {
73 ss << dims[1];
74 if (num > 2) {
75 ss << "," << dims[2];
76 if (num > 3) {
77 ss << "," << dims[3];
78 }
79 }
80 }
81 ss << ").";
82 const std::string errorMsg(ss.str());
83 throw escript::ValueError(errorMsg);
84 }
85
86 void Assemble_PDE(const NodeFile* nodes, const ElementFile* elements,
87 escript::ASM_ptr S, escript::Data& F,
88 const escript::Data& A, const escript::Data& B,
89 const escript::Data& C, const escript::Data& D,
90 const escript::Data& X, const escript::Data& Y)
91 {
92 if (!nodes || !elements || (S==NULL && F.isEmpty()))
93 return;
94
95 if (F.isEmpty() && (!X.isEmpty() || !Y.isEmpty())) {
96 throw escript::ValueError("Assemble_PDE: right hand side coefficients are non-zero but no right hand side vector given.");
97 }
98
99 if (S==NULL && !A.isEmpty() && !B.isEmpty() && !C.isEmpty() && !D.isEmpty()) {
100 throw escript::ValueError("Assemble_PDE: coefficients are non-zero but no matrix is given.");
101 }
102
103 // get the function space for this assemblage call
104 int funcspace = -1;
105 if (!A.isEmpty()) funcspace=A.getFunctionSpace().getTypeCode();
106 if (!B.isEmpty()) funcspace=B.getFunctionSpace().getTypeCode();
107 if (!C.isEmpty()) funcspace=C.getFunctionSpace().getTypeCode();
108 if (!D.isEmpty()) funcspace=D.getFunctionSpace().getTypeCode();
109 if (!X.isEmpty()) funcspace=X.getFunctionSpace().getTypeCode();
110 if (!Y.isEmpty()) funcspace=Y.getFunctionSpace().getTypeCode();
111 if (funcspace == -1)
112 return; // all data are empty
113
114 // check if all function spaces are the same
115 if (!A.isEmpty() && A.getFunctionSpace().getTypeCode()!=funcspace) {
116 throw escript::ValueError("Assemble_PDE: unexpected function space type for coefficient A");
117 } else if (!B.isEmpty() && B.getFunctionSpace().getTypeCode()!=funcspace) {
118 throw escript::ValueError("Assemble_PDE: unexpected function space type for coefficient B");
119 } else if (!C.isEmpty() && C.getFunctionSpace().getTypeCode()!=funcspace) {
120 throw escript::ValueError("Assemble_PDE: unexpected function space type for coefficient C");
121 } else if (!D.isEmpty() && D.getFunctionSpace().getTypeCode()!=funcspace) {
122 throw escript::ValueError("Assemble_PDE: unexpected function space type for coefficient D");
123 } else if (!X.isEmpty() && X.getFunctionSpace().getTypeCode()!=funcspace) {
124 throw escript::ValueError("Assemble_PDE: unexpected function space type for coefficient X");
125 } else if (!Y.isEmpty() && Y.getFunctionSpace().getTypeCode()!=funcspace) {
126 throw escript::ValueError("Assemble_PDE: unexpected function space type for coefficient Y");
127 }
128
129 // get value type
130 bool isComplex = false;
131 isComplex = isComplex || (!A.isEmpty() && A.isComplex());
132 isComplex = isComplex || (!B.isEmpty() && B.isComplex());
133 isComplex = isComplex || (!C.isEmpty() && C.isComplex());
134 isComplex = isComplex || (!D.isEmpty() && D.isComplex());
135 isComplex = isComplex || (!X.isEmpty() && X.isComplex());
136 isComplex = isComplex || (!Y.isEmpty() && Y.isComplex());
137
138 bool reducedIntegrationOrder;
139 if (funcspace==FINLEY_ELEMENTS) {
140 reducedIntegrationOrder=false;
141 } else if (funcspace==FINLEY_FACE_ELEMENTS) {
142 reducedIntegrationOrder=false;
143 } else if (funcspace==FINLEY_CONTACT_ELEMENTS_1) {
144 reducedIntegrationOrder=false;
145 } else if (funcspace==FINLEY_CONTACT_ELEMENTS_2) {
146 reducedIntegrationOrder=false;
147 } else if (funcspace==FINLEY_REDUCED_ELEMENTS) {
148 reducedIntegrationOrder=true;
149 } else if (funcspace==FINLEY_REDUCED_FACE_ELEMENTS) {
150 reducedIntegrationOrder=true;
151 } else if (funcspace==FINLEY_REDUCED_CONTACT_ELEMENTS_1) {
152 reducedIntegrationOrder=true;
153 } else if (funcspace==FINLEY_REDUCED_CONTACT_ELEMENTS_2) {
154 reducedIntegrationOrder=true;
155 } else if (funcspace==FINLEY_POINTS) {
156 reducedIntegrationOrder=false;
157 } else {
158 throw escript::ValueError("Assemble_PDE: assemblage failed because of illegal function space.");
159 return;
160 }
161
162 // get assemblage parameters
163 AssembleParameters p(nodes, elements, S, F, reducedIntegrationOrder);
164
165 // check if sample numbers are the same
166 if (!A.numSamplesEqual(p.numQuadTotal, elements->numElements)) {
167 setNumSamplesError("A", p.numQuadTotal, elements->numElements);
168 } else if (!B.numSamplesEqual(p.numQuadTotal, elements->numElements)) {
169 setNumSamplesError("B", p.numQuadTotal, elements->numElements);
170 } else if (!C.numSamplesEqual(p.numQuadTotal, elements->numElements)) {
171 setNumSamplesError("C", p.numQuadTotal, elements->numElements);
172 } else if (!D.numSamplesEqual(p.numQuadTotal, elements->numElements)) {
173 setNumSamplesError("D", p.numQuadTotal, elements->numElements);
174 } else if (!X.numSamplesEqual(p.numQuadTotal, elements->numElements)) {
175 setNumSamplesError("X", p.numQuadTotal, elements->numElements);
176 } else if (!Y.numSamplesEqual(p.numQuadTotal, elements->numElements)) {
177 setNumSamplesError("Y", p.numQuadTotal, elements->numElements);
178 }
179
180 // check the dimensions:
181 if (p.numEqu != p. numComp) {
182 throw escript::ValueError("Assemble_PDE requires number of equations == number of solutions.");
183 } else if (p.numEqu == 1) {
184 const int dimensions[2] = { p.numDim, p.numDim };
185 if (!A.isDataPointShapeEqual(2, dimensions)) {
186 setShapeError("A", 2, dimensions);
187 } else if (!B.isDataPointShapeEqual(1, dimensions)) {
188 setShapeError("B", 1, dimensions);
189 } else if (!C.isDataPointShapeEqual(1, dimensions)) {
190 setShapeError("C", 1, dimensions);
191 } else if (!D.isDataPointShapeEqual(0, dimensions)) {
192 throw escript::ValueError("Assemble_PDE: coefficient D must be rank 0.");
193 } else if (!X.isDataPointShapeEqual(1, dimensions)) {
194 setShapeError("X", 1, dimensions);
195 } else if (!Y.isDataPointShapeEqual(0, dimensions)) {
196 throw escript::ValueError("Assemble_PDE: coefficient Y must be rank 0.");
197 }
198 } else {
199 const int dimAB[4] = { p.numEqu, p.numDim, p.numComp, p.numDim };
200 const int dimCD[3] = { p.numEqu, p.numComp, p.numDim };
201 if (!A.isDataPointShapeEqual(4, dimAB)) {
202 setShapeError("A", 4, dimAB);
203 } else if (!B.isDataPointShapeEqual(3, dimAB)) {
204 setShapeError("B", 3, dimAB);
205 } else if (!C.isDataPointShapeEqual(3, dimCD)) {
206 setShapeError("C", 3, dimCD);
207 } else if (!D.isDataPointShapeEqual(2, dimCD)) {
208 setShapeError("D", 2, dimCD);
209 } else if (!X.isDataPointShapeEqual(2, dimAB)) {
210 setShapeError("X", 2, dimAB);
211 } else if (!Y.isDataPointShapeEqual(1, dimAB)) {
212 setShapeError("Y", 1, dimAB);
213 }
214 }
215
216 if (p.numSides == 1) {
217 if (funcspace == FINLEY_POINTS) {
218 if (!A.isEmpty() || !B.isEmpty() || !C.isEmpty() || !X.isEmpty()) {
219 throw escript::ValueError("Assemble_PDE: Point elements require A, B, C and X to be empty.");
220 } else {
221 if (isComplex) {
222 Assemble_PDE_Points<cplx_t>(p, D, Y);
223 } else {
224 Assemble_PDE_Points<real_t>(p, D, Y);
225 }
226 }
227 } else if (p.numEqu > 1) { // system of PDEs
228 if (p.numDim == 3) {
229 if (isComplex) {
230 Assemble_PDE_System_3D<cplx_t>(p, A, B, C, D, X, Y);
231 } else {
232 Assemble_PDE_System_3D<real_t>(p, A, B, C, D, X, Y);
233 }
234 } else if (p.numDim == 2) {
235 if (isComplex) {
236 Assemble_PDE_System_2D<cplx_t>(p, A, B, C, D, X, Y);
237 } else {
238 Assemble_PDE_System_2D<real_t>(p, A, B, C, D, X, Y);
239 }
240 } else if (p.numDim == 1) {
241 Assemble_PDE_System_1D(p, A, B, C, D, X, Y);
242 } else {
243 throw escript::ValueError("Assemble_PDE supports spatial dimensions 1,2,3 only.");
244 }
245 } else { // single PDE
246 if (p.numDim == 3) {
247 if (isComplex) {
248 Assemble_PDE_Single_3D<cplx_t>(p, A, B, C, D, X, Y);
249 } else {
250 Assemble_PDE_Single_3D<real_t>(p, A, B, C, D, X, Y);
251 }
252 } else if (p.numDim == 2) {
253 if (isComplex) {
254 Assemble_PDE_Single_2D<cplx_t>(p, A, B, C, D, X, Y);
255 } else {
256 Assemble_PDE_Single_2D<real_t>(p, A, B, C, D, X, Y);
257 }
258 } else if (p.numDim == 1) {
259 Assemble_PDE_Single_1D(p, A, B, C, D, X, Y);
260 } else {
261 throw escript::ValueError("Assemble_PDE supports spatial dimensions 1,2,3 only.");
262 }
263 }
264 } else if (p.numSides == 2) {
265 if (!A.isEmpty() || !B.isEmpty() || !C.isEmpty() || !X.isEmpty()) {
266 throw escript::ValueError("Assemble_PDE: Contact elements require A, B, C and X to be empty.");
267 } else if (p.numEqu > 1) { // system of PDEs
268 if (isComplex) {
269 Assemble_PDE_System_C<cplx_t>(p, D, Y);
270 } else {
271 Assemble_PDE_System_C<real_t>(p, D, Y);
272 }
273 } else { // single PDE
274 if (isComplex) {
275 Assemble_PDE_Single_C<cplx_t>(p, D, Y);
276 } else {
277 Assemble_PDE_Single_C<real_t>(p, D, Y);
278 }
279 }
280 } else {
281 throw escript::ValueError("Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
282 }
283 }
284
285 } // namespace finley
286

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26