/[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 6009 - (show annotations)
Wed Mar 2 04:13:26 2016 UTC (3 years, 1 month ago) by caltinay
File size: 17027 byte(s)
Much needed sync with trunk...

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

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