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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6079 - (show annotations)
Mon Mar 21 12:22:38 2016 UTC (2 years, 10 months ago) by caltinay
File size: 16645 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 #include "Assemble.h"
18 #include "ShapeTable.h"
19 #include "Util.h"
20
21 namespace dudley {
22
23 void Assemble_LumpedSystem(const NodeFile* nodes, const ElementFile* elements,
24 escript::Data& lumpedMat, const escript::Data& D,
25 bool useHRZ)
26 {
27 if (!nodes || !elements || lumpedMat.isEmpty() || D.isEmpty())
28 return;
29
30 const int funcspace = D.getFunctionSpace().getTypeCode();
31 bool reducedIntegrationOrder;
32 // check function space of D
33 if (funcspace == DUDLEY_ELEMENTS) {
34 reducedIntegrationOrder = false;
35 } else if (funcspace == DUDLEY_FACE_ELEMENTS) {
36 reducedIntegrationOrder = false;
37 } else if (funcspace == DUDLEY_REDUCED_ELEMENTS) {
38 reducedIntegrationOrder = true;
39 } else if (funcspace == DUDLEY_REDUCED_FACE_ELEMENTS) {
40 reducedIntegrationOrder = true;
41 } else {
42 throw DudleyException("Assemble_LumpedSystem: assemblage failed because of illegal function space.");
43 }
44
45 // initialize parameters
46 AssembleParameters p;
47 Assemble_getAssembleParameters(nodes, elements, escript::ASM_ptr(),
48 lumpedMat, reducedIntegrationOrder, &p);
49
50 // check if all function spaces are the same
51 if (!D.numSamplesEqual(p.numQuad, elements->numElements)) {
52 std::stringstream ss;
53 ss << "Assemble_LumpedSystem: sample points of coefficient D "
54 "don't match (" << p.numQuad << ","
55 << elements->numElements << ")";
56 const std::string msg(ss.str());
57 throw DudleyException(msg);
58 }
59
60 // check the dimensions
61 if (p.numEqu == 1) {
62 const escript::DataTypes::ShapeType dimensions; //dummy
63 if (D.getDataPointShape() != dimensions) {
64 throw DudleyException("Assemble_LumpedSystem: coefficient D, rank 0 expected.");
65 }
66 } else {
67 const escript::DataTypes::ShapeType dimensions(1, p.numEqu);
68 if (D.getDataPointShape() != dimensions) {
69 std::stringstream ss;
70 ss << "Assemble_LumpedSystem: coefficient D, expected "
71 "shape (" << p.numEqu << ",)";
72 const std::string msg(ss.str());
73 throw DudleyException(msg);
74 }
75 }
76
77 lumpedMat.requireWrite();
78 double* lumpedMat_p = lumpedMat.getSampleDataRW(0);
79
80 if (funcspace==DUDLEY_POINTS) {
81 #pragma omp parallel
82 {
83 for (int color=elements->minColor; color<=elements->maxColor; color++) {
84 // loop over all elements
85 #pragma omp for
86 for (index_t e=0; e<elements->numElements; e++) {
87 if (elements->Color[e]==color) {
88 const double* D_p = D.getSampleDataRO(e);
89 util::addScatter(1,
90 &p.row_DOF[elements->Nodes[INDEX2(0,e,p.NN)]],
91 p.numEqu, D_p, lumpedMat_p,
92 p.row_DOF_UpperBound);
93 } // end color check
94 } // end element loop
95 } // end color loop
96 } // end parallel region
97 } else {
98 bool expandedD = D.actsExpanded();
99 const double *S = NULL;
100 if (!getQuadShape(elements->numDim, reducedIntegrationOrder, &S)) {
101 throw DudleyException("Assemble_LumpedSystem: Unable to locate shape function.");
102 }
103 #pragma omp parallel
104 {
105 std::vector<double> EM_lumpedMat(p.numShapes * p.numEqu);
106 std::vector<index_t> row_index(p.numShapes);
107
108 if (p.numEqu == 1) { // single equation
109 if (expandedD) { // with expanded D
110 for (int color = elements->minColor; color <= elements->maxColor; color++) {
111 // loop over all elements
112 #pragma omp for
113 for (index_t e = 0; e < elements->numElements; e++) {
114 if (elements->Color[e] == color) {
115 const double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
116 const double* D_p = D.getSampleDataRO(e);
117 if (useHRZ) {
118 double m_t = 0; // mass of the element
119 for (int q = 0; q < p.numQuad; q++)
120 m_t += vol * D_p[INDEX2(q, 0, p.numQuad)];
121 double diagS = 0; // diagonal sum
122 double rtmp;
123 for (int s = 0; s < p.numShapes; s++) {
124 rtmp = 0.;
125 for (int q = 0; q < p.numQuad; q++)
126 rtmp +=
127 vol * D_p[INDEX2(q, 0, p.numQuad)] * S[INDEX2(s, q, p.numShapes)] *
128 S[INDEX2(s, q, p.numShapes)];
129 EM_lumpedMat[INDEX2(0, s, p.numEqu)] = rtmp;
130 diagS += rtmp;
131 }
132 // rescale diagonals by m_t/diagS to ensure
133 // consistent mass over element
134 rtmp = m_t / diagS;
135 for (int s = 0; s < p.numShapes; s++)
136 EM_lumpedMat[INDEX2(0, s, p.numEqu)] *= rtmp;
137 } else { // row-sum lumping
138 for (int s = 0; s < p.numShapes; s++) {
139 double rtmp = 0.;
140 for (int q = 0; q < p.numQuad; q++)
141 rtmp += vol * S[INDEX2(s, q, p.numShapes)] * D_p[INDEX2(q, 0, p.numQuad)];
142 EM_lumpedMat[INDEX2(0, s, p.numEqu)] = rtmp;
143 }
144 }
145 for (int q = 0; q < p.numShapes; q++)
146 row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
147 util::addScatter(p.numShapes, &row_index[0],
148 p.numEqu, &EM_lumpedMat[0], lumpedMat_p,
149 p.row_DOF_UpperBound);
150 } // end color check
151 } // end element loop
152 } // end color loop
153 } else { // with constant D
154 for (int color = elements->minColor; color <= elements->maxColor; color++) {
155 // loop over all elements
156 #pragma omp for
157 for (index_t e = 0; e < elements->numElements; e++) {
158 if (elements->Color[e] == color) {
159 const double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
160 const double* D_p = D.getSampleDataRO(e);
161 if (useHRZ) { // HRZ lumping
162 // mass of the element
163 const double m_t = vol*p.numQuad;
164 double diagS = 0; // diagonal sum
165 double rtmp;
166 for (int s = 0; s < p.numShapes; s++) {
167 rtmp = 0.;
168 for (int q = 0; q < p.numQuad; q++) {
169 rtmp += vol * S[INDEX2(s, q, p.numShapes)] * S[INDEX2(s, q, p.numShapes)];
170 }
171 EM_lumpedMat[INDEX2(0, s, p.numEqu)] = rtmp;
172 diagS += rtmp;
173 }
174 // rescale diagonals by m_t/diagS to ensure
175 // consistent mass over element
176 rtmp = m_t / diagS * D_p[0];
177 for (int s = 0; s < p.numShapes; s++)
178 EM_lumpedMat[INDEX2(0, s, p.numEqu)] *= rtmp;
179 } else { // row-sum lumping
180 for (int s = 0; s < p.numShapes; s++) {
181 double rtmp = 0.;
182 for (int q = 0; q < p.numQuad; q++)
183 rtmp += vol * S[INDEX2(s, q, p.numShapes)];
184 EM_lumpedMat[INDEX2(0, s, p.numEqu)] = rtmp * D_p[0];
185 }
186 }
187 for (int q = 0; q < p.numShapes; q++)
188 row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
189 util::addScatter(p.numShapes, &row_index[0],
190 p.numEqu, &EM_lumpedMat[0], lumpedMat_p,
191 p.row_DOF_UpperBound);
192 } // end color check
193 } // end element loop
194 } // end color loop
195 }
196
197 } else { // system of equations
198 if (expandedD) { // with expanded D
199 for (int color = elements->minColor; color <= elements->maxColor; color++) {
200 // loop over all elements
201 #pragma omp for
202 for (index_t e = 0; e < elements->numElements; e++) {
203 if (elements->Color[e] == color) {
204 const double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
205 const double* D_p = D.getSampleDataRO(e);
206
207 if (useHRZ) { // HRZ lumping
208 for (int k = 0; k < p.numEqu; k++) {
209 double m_t = 0; // mass of the element
210 for (int q = 0; q < p.numQuad; q++)
211 m_t += vol * D_p[INDEX3(k, q, 0, p.numEqu, p.numQuad)];
212
213 double diagS = 0; // diagonal sum
214 double rtmp;
215 for (int s = 0; s < p.numShapes; s++) {
216 rtmp = 0.;
217 for (int q = 0; q < p.numQuad; q++)
218 rtmp +=
219 vol * D_p[INDEX3(k, q, 0, p.numEqu, p.numQuad)] *
220 S[INDEX2(s, q, p.numShapes)] * S[INDEX2(s, q, p.numShapes)];
221 EM_lumpedMat[INDEX2(k, s, p.numEqu)] = rtmp;
222 diagS += rtmp;
223 }
224 // rescale diagonals by m_t/diagS to
225 // ensure consistent mass over element
226 rtmp = m_t / diagS;
227 for (int s = 0; s < p.numShapes; s++)
228 EM_lumpedMat[INDEX2(k, s, p.numEqu)] *= rtmp;
229 }
230 } else { // row-sum lumping
231 for (int s = 0; s < p.numShapes; s++) {
232 for (int k = 0; k < p.numEqu; k++) {
233 double rtmp = 0.;
234 for (int q = 0; q < p.numQuad; q++)
235 rtmp +=
236 vol * S[INDEX2(s, q, p.numShapes)] *
237 D_p[INDEX3(k, q, 0, p.numEqu, p.numQuad)];
238 EM_lumpedMat[INDEX2(k, s, p.numEqu)] = rtmp;
239 }
240 }
241 }
242 for (int q = 0; q < p.numShapes; q++)
243 row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
244 util::addScatter(p.numShapes, &row_index[0],
245 p.numEqu, &EM_lumpedMat[0], lumpedMat_p,
246 p.row_DOF_UpperBound);
247 } // end color check
248 } // end element loop
249 } // end color loop
250 } else { // with constant D
251 for (int color = elements->minColor; color <= elements->maxColor; color++) {
252 // loop over all elements
253 #pragma omp for
254 for (index_t e = 0; e < elements->numElements; e++) {
255 if (elements->Color[e] == color) {
256 const double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
257 const double* D_p = D.getSampleDataRO(e);
258
259 if (useHRZ) { // HRZ lumping
260 double m_t = vol * p.numQuad; // mass of the element
261 double diagS = 0; // diagonal sum
262 double rtmp;
263 for (int s = 0; s < p.numShapes; s++) {
264 rtmp = 0.;
265 for (int q = 0; q < p.numQuad; q++)
266 rtmp += vol * S[INDEX2(s, q, p.numShapes)] * S[INDEX2(s, q, p.numShapes)];
267 for (int k = 0; k < p.numEqu; k++)
268 EM_lumpedMat[INDEX2(k, s, p.numEqu)] = rtmp;
269 diagS += rtmp;
270 }
271 // rescale diagonals by m_t/diagS to ensure
272 // consistent mass over element
273 rtmp = m_t / diagS;
274 for (int s = 0; s < p.numShapes; s++) {
275 for (int k = 0; k < p.numEqu; k++)
276 EM_lumpedMat[INDEX2(k, s, p.numEqu)] *= rtmp * D_p[k];
277 }
278 } else { // row-sum lumping
279 for (int s = 0; s < p.numShapes; s++) {
280 for (int k = 0; k < p.numEqu; k++) {
281 double rtmp = 0.;
282 for (int q = 0; q < p.numQuad; q++)
283 rtmp += vol * S[INDEX2(s, q, p.numShapes)];
284 EM_lumpedMat[INDEX2(k, s, p.numEqu)] = rtmp * D_p[k];
285 }
286 }
287 }
288 for (int q = 0; q < p.numShapes; q++)
289 row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
290 util::addScatter(p.numShapes, &row_index[0],
291 p.numEqu, &EM_lumpedMat[0], lumpedMat_p,
292 p.row_DOF_UpperBound);
293 } // end color check
294 } // end element loop
295 } // end color loop
296 }
297 }
298 } // end parallel region
299 }
300 }
301
302 } // namespace dudley
303

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26