/[escript]/branches/trilinos_from_5897/dudley/src/Assemble_PDE.cpp
ViewVC logotype

Contents of /branches/trilinos_from_5897/dudley/src/Assemble_PDE.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6090 - (show annotations)
Wed Mar 23 06:35:54 2016 UTC (3 years, 1 month ago) by caltinay
File size: 9150 byte(s)
Simplified dudley PDE routines.

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2016 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
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 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17 /****************************************************************************
18
19 Assembles the system of numEqu PDEs into the stiffness matrix S and right
20 hand side F:
21
22 -div(A*grad u)-div(B*u)+C*grad u + D*u = -div X + Y
23
24 -(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
25
26 u has numEqu components.
27 Shape of the coefficients:
28
29 A = numEqu x numDim x numEqu x numDim
30 B = numDim x numEqu x numEqu
31 C = numEqu x numDim x numEqu
32 D = numEqu x numEqu
33 X = numEqu x numDim
34 Y = numEqu
35
36 The coefficients A,B,C,D,X and Y have to be defined on the integration points
37 or not present (i.e. empty).
38
39 S and F have to be initialised before the routine is called. S or F can be
40 NULL. In this case the left or the right hand side of the PDE is not
41 processed.
42
43 The routine does not consider any boundary conditions.
44
45 *****************************************************************************/
46
47 #include "Assemble.h"
48 #include "Util.h"
49
50 namespace dudley {
51
52 inline void setNumSamplesError(const char* c, int n0, int n1)
53 {
54 std::stringstream ss;
55 ss << "Assemble_PDE: number of sample points of coefficient " << c
56 << " don't match (" << n0 << "," << n1 << ").";
57 const std::string errorMsg(ss.str());
58 throw DudleyException(errorMsg);
59 }
60
61 inline void setShapeError(const char* c, int num, const int* dims)
62 {
63 std::stringstream ss;
64 ss << "Assemble_PDE: shape of coefficient " << c
65 << " does not match (" << dims[0] << ",";
66 if (num > 1) {
67 ss << dims[1];
68 if (num > 2) {
69 ss << "," << dims[2];
70 if (num > 3) {
71 ss << "," << dims[3];
72 }
73 }
74 }
75 ss << ").";
76 const std::string errorMsg(ss.str());
77 throw DudleyException(errorMsg);
78 }
79
80 void Assemble_PDE(const NodeFile* nodes, const ElementFile* elements,
81 escript::ASM_ptr S, escript::Data& F,
82 const escript::Data& A, const escript::Data& B,
83 const escript::Data& C, const escript::Data& D,
84 const escript::Data& X, const escript::Data& Y)
85 {
86 if (!nodes || !elements || (S.get()==NULL && F.isEmpty()))
87 return;
88
89 if (F.isEmpty() && (!X.isEmpty() || !Y.isEmpty())) {
90 throw DudleyException("Assemble_PDE: right hand side coefficients are non-zero but no right hand side vector given.");
91 }
92
93 if (S.get()==NULL && !A.isEmpty() && !B.isEmpty() && !C.isEmpty() && !D.isEmpty()) {
94 throw DudleyException("Assemble_PDE: coefficients are non-zero but no matrix is given.");
95 }
96
97 // get the functionspace for this assemblage call
98 int funcspace = -1;
99 if (!A.isEmpty()) funcspace=A.getFunctionSpace().getTypeCode();
100 if (!B.isEmpty()) funcspace=B.getFunctionSpace().getTypeCode();
101 if (!C.isEmpty()) funcspace=C.getFunctionSpace().getTypeCode();
102 if (!D.isEmpty()) funcspace=D.getFunctionSpace().getTypeCode();
103 if (!X.isEmpty()) funcspace=X.getFunctionSpace().getTypeCode();
104 if (!Y.isEmpty()) funcspace=Y.getFunctionSpace().getTypeCode();
105 if (funcspace == -1)
106 return; // all data are empty
107
108 // check if all function spaces are the same
109 if (!A.isEmpty() && A.getFunctionSpace().getTypeCode()!=funcspace) {
110 throw DudleyException("Assemble_PDE: unexpected function space type for coefficient A");
111 } else if (!B.isEmpty() && B.getFunctionSpace().getTypeCode()!=funcspace) {
112 throw DudleyException("Assemble_PDE: unexpected function space type for coefficient B");
113 } else if (!C.isEmpty() && C.getFunctionSpace().getTypeCode()!=funcspace) {
114 throw DudleyException("Assemble_PDE: unexpected function space type for coefficient C");
115 } else if (!D.isEmpty() && D.getFunctionSpace().getTypeCode()!=funcspace) {
116 throw DudleyException("Assemble_PDE: unexpected function space type for coefficient D");
117 } else if (!X.isEmpty() && X.getFunctionSpace().getTypeCode()!=funcspace) {
118 throw DudleyException("Assemble_PDE: unexpected function space type for coefficient X");
119 } else if (!Y.isEmpty() && Y.getFunctionSpace().getTypeCode()!=funcspace) {
120 throw DudleyException("Assemble_PDE: unexpected function space type for coefficient Y");
121 }
122
123 bool reducedIntegrationOrder;
124 if (funcspace == DUDLEY_ELEMENTS) {
125 reducedIntegrationOrder = false;
126 } else if (funcspace == DUDLEY_FACE_ELEMENTS) {
127 reducedIntegrationOrder = false;
128 } else if (funcspace == DUDLEY_REDUCED_ELEMENTS) {
129 reducedIntegrationOrder = true;
130 } else if (funcspace == DUDLEY_REDUCED_FACE_ELEMENTS) {
131 reducedIntegrationOrder = true;
132 } else if (funcspace == DUDLEY_POINTS) {
133 reducedIntegrationOrder = true;
134 } else {
135 throw DudleyException("Assemble_PDE: assemblage failed because of illegal function space.");
136 }
137
138 // get assemblage parameters
139 AssembleParameters p(nodes, elements, S, F, reducedIntegrationOrder);
140
141 // check if sample numbers are the same
142 if (!A.numSamplesEqual(p.numQuad, elements->numElements)) {
143 setNumSamplesError("A", p.numQuad, elements->numElements);
144 } else if (!B.numSamplesEqual(p.numQuad, elements->numElements)) {
145 setNumSamplesError("B", p.numQuad, elements->numElements);
146 } else if (!C.numSamplesEqual(p.numQuad, elements->numElements)) {
147 setNumSamplesError("C", p.numQuad, elements->numElements);
148 } else if (!D.numSamplesEqual(p.numQuad, elements->numElements)) {
149 setNumSamplesError("D", p.numQuad, elements->numElements);
150 } else if (!X.numSamplesEqual(p.numQuad, elements->numElements)) {
151 setNumSamplesError("X", p.numQuad, elements->numElements);
152 } else if (!Y.numSamplesEqual(p.numQuad, elements->numElements)) {
153 setNumSamplesError("Y", p.numQuad, elements->numElements);
154 }
155
156 // check the dimensions
157 if (p.numEqu == 1) {
158 const int dimensions[2] = { p.numDim, p.numDim };
159 if (!A.isDataPointShapeEqual(2, dimensions)) {
160 setShapeError("A", 2, dimensions);
161 } else if (!B.isDataPointShapeEqual(1, dimensions)) {
162 setShapeError("B", 1, dimensions);
163 } else if (!C.isDataPointShapeEqual(1, dimensions)) {
164 setShapeError("C", 1, dimensions);
165 } else if (!D.isDataPointShapeEqual(0, dimensions)) {
166 throw DudleyException("Assemble_PDE: coefficient D must be rank 0.");
167 } else if (!X.isDataPointShapeEqual(1, dimensions)) {
168 setShapeError("X", 1, dimensions);
169 } else if (!Y.isDataPointShapeEqual(0, dimensions)) {
170 throw DudleyException("Assemble_PDE: coefficient Y must be rank 0.");
171 }
172 } else {
173 const int dimAB[4] = { p.numEqu, p.numDim, p.numEqu, p.numDim };
174 const int dimCD[3] = { p.numEqu, p.numEqu, p.numDim };
175 if (!A.isDataPointShapeEqual(4, dimAB)) {
176 setShapeError("A", 4, dimAB);
177 } else if (!B.isDataPointShapeEqual(3, dimAB)) {
178 setShapeError("B", 3, dimAB);
179 } else if (!C.isDataPointShapeEqual(3, dimCD)) {
180 setShapeError("C", 3, dimCD);
181 } else if (!D.isDataPointShapeEqual(2, dimCD)) {
182 setShapeError("D", 2, dimCD);
183 } else if (!X.isDataPointShapeEqual(2, dimAB)) {
184 setShapeError("X", 2, dimAB);
185 } else if (!Y.isDataPointShapeEqual(1, dimAB)) {
186 setShapeError("Y", 1, dimAB);
187 }
188 }
189
190 if (funcspace==DUDLEY_POINTS) {
191 if (!A.isEmpty() || !B.isEmpty() || !C.isEmpty() || !X.isEmpty()) {
192 throw DudleyException("Dudley_Assemble_PDE: Point elements require A, B, C and X to be empty.");
193 } else {
194 Assemble_PDE_Points(p, D, Y);
195 }
196 } else {
197 if (p.numEqu > 1) {
198 // system of PDEs
199 if (p.numDim == 3) {
200 Assemble_PDE_System_3D(p, A, B, C, D, X, Y);
201 } else if (p.numDim == 2) {
202 Assemble_PDE_System_2D(p, A, B, C, D, X, Y);
203 } else {
204 throw DudleyException("Assemble_PDE supports spatial dimensions 2 and 3 only.");
205 }
206 } else {
207 // single PDE
208 if (p.numDim == 3) {
209 Assemble_PDE_Single_3D(p, A, B, C, D, X, Y);
210 } else if (p.numDim == 2) {
211 Assemble_PDE_Single_2D(p, A, B, C, D, X, Y);
212 } else {
213 throw DudleyException("Assemble_PDE supports spatial dimensions 2 and 3 only.");
214 }
215 }
216 }
217 }
218
219 } // namespace dudley
220

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26