/[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 6079 - (show annotations)
Mon Mar 21 12:22:38 2016 UTC (2 years, 11 months ago) by caltinay
File size: 9417 byte(s)
Big commit - making dudley much more like finley to make it more
managable. Fixed quite a few issues that had been fixed in finley.
Disposed of all ReducedNode/ReducedDOF entities that dudley never supported.
Compiles and passes tests.

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 numComp components.
27 Shape of the coefficients:
28
29 A = numEqu x numDim x numComp x numDim
30 B = numDim x numEqu x numComp
31 C = numEqu x numDim x numComp
32 D = numEqu x numComp
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;
140 Assemble_getAssembleParameters(nodes, elements, S, F, reducedIntegrationOrder, &p);
141
142 // check if sample numbers are the same
143 if (!A.numSamplesEqual(p.numQuad, elements->numElements)) {
144 setNumSamplesError("A", p.numQuad, elements->numElements);
145 } else if (!B.numSamplesEqual(p.numQuad, elements->numElements)) {
146 setNumSamplesError("B", p.numQuad, elements->numElements);
147 } else if (!C.numSamplesEqual(p.numQuad, elements->numElements)) {
148 setNumSamplesError("C", p.numQuad, elements->numElements);
149 } else if (!D.numSamplesEqual(p.numQuad, elements->numElements)) {
150 setNumSamplesError("D", p.numQuad, elements->numElements);
151 } else if (!X.numSamplesEqual(p.numQuad, elements->numElements)) {
152 setNumSamplesError("X", p.numQuad, elements->numElements);
153 } else if (!Y.numSamplesEqual(p.numQuad, elements->numElements)) {
154 setNumSamplesError("Y", p.numQuad, elements->numElements);
155 }
156
157 // check the dimensions
158 if (p.numEqu != p.numComp) {
159 throw DudleyException("Assemble_PDE requires number of equations == number of solutions.");
160 } else if (p.numEqu == 1) {
161 const int dimensions[2] = { p.numDim, p.numDim };
162 if (!A.isDataPointShapeEqual(2, dimensions)) {
163 setShapeError("A", 2, dimensions);
164 } else if (!B.isDataPointShapeEqual(1, dimensions)) {
165 setShapeError("B", 1, dimensions);
166 } else if (!C.isDataPointShapeEqual(1, dimensions)) {
167 setShapeError("C", 1, dimensions);
168 } else if (!D.isDataPointShapeEqual(0, dimensions)) {
169 throw DudleyException("Assemble_PDE: coefficient D must be rank 0.");
170 } else if (!X.isDataPointShapeEqual(1, dimensions)) {
171 setShapeError("X", 1, dimensions);
172 } else if (!Y.isDataPointShapeEqual(0, dimensions)) {
173 throw DudleyException("Assemble_PDE: coefficient Y must be rank 0.");
174 }
175 } else {
176 const int dimAB[4] = { p.numEqu, p.numDim, p.numComp, p.numDim };
177 const int dimCD[3] = { p.numEqu, p.numComp, p.numDim };
178 if (!A.isDataPointShapeEqual(4, dimAB)) {
179 setShapeError("A", 4, dimAB);
180 } else if (!B.isDataPointShapeEqual(3, dimAB)) {
181 setShapeError("B", 3, dimAB);
182 } else if (!C.isDataPointShapeEqual(3, dimCD)) {
183 setShapeError("C", 3, dimCD);
184 } else if (!D.isDataPointShapeEqual(2, dimCD)) {
185 setShapeError("D", 2, dimCD);
186 } else if (!X.isDataPointShapeEqual(2, dimAB)) {
187 setShapeError("X", 2, dimAB);
188 } else if (!Y.isDataPointShapeEqual(1, dimAB)) {
189 setShapeError("Y", 1, dimAB);
190 }
191 }
192
193 if (funcspace==DUDLEY_POINTS) {
194 if (!A.isEmpty() || !B.isEmpty() || !C.isEmpty() || !X.isEmpty()) {
195 throw DudleyException("Dudley_Assemble_PDE: Point elements require A, B, C and X to be empty.");
196 } else {
197 Assemble_PDE_Points(p, elements, S, F, D, Y);
198 }
199 } else {
200 if (p.numEqu > 1) {
201 // system of PDEs
202 if (p.numDim == 3) {
203 Assemble_PDE_System_3D(p, elements, S, F, A, B, C, D, X, Y);
204 } else if (p.numDim == 2) {
205 Assemble_PDE_System_2D(p, elements, S, F, A, B, C, D, X, Y);
206 } else {
207 throw DudleyException("Assemble_PDE supports spatial dimensions 2 and 3 only.");
208 }
209 } else {
210 // single PDE
211 if (p.numDim == 3) {
212 Assemble_PDE_Single_3D(p, elements, S, F, A, B, C, D, X, Y);
213 } else if (p.numDim == 2) {
214 Assemble_PDE_Single_2D(p, elements, S, F, A, B, C, D, X, Y);
215 } else {
216 throw DudleyException("Assemble_PDE supports spatial dimensions 2 and 3 only.");
217 }
218 }
219 }
220 }
221
222 } // namespace dudley
223

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