/[escript]/trunk/finley/src/Assemble_LumpedSystem.c
ViewVC logotype

Contents of /trunk/finley/src/Assemble_LumpedSystem.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3379 - (show annotations)
Wed Nov 24 04:48:49 2010 UTC (8 years, 1 month ago) by gross
File MIME type: text/plain
File size: 15065 byte(s)
some clarification on lumping
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14
15 /**************************************************************/
16
17 /* assembles the mass matrix in lumped form */
18
19 /* The coefficient D has to be defined on the integration points or not present. */
20
21 /* lumpedMat has to be initialized before the routine is called. */
22
23 /**************************************************************/
24
25 #include "Assemble.h"
26 #include "Util.h"
27 #ifdef _OPENMP
28 #include <omp.h>
29 #endif
30
31
32 /**************************************************************/
33
34 void Finley_Assemble_LumpedSystem(Finley_NodeFile* nodes,Finley_ElementFile* elements, escriptDataC* lumpedMat, escriptDataC* D, const bool_t useHRZ)
35 {
36
37 bool_t reducedIntegrationOrder=FALSE, expandedD;
38 char error_msg[LenErrorMsg_MAX];
39 Finley_Assemble_Parameters p;
40 dim_t dimensions[ESCRIPT_MAX_DATA_RANK], k, e, len_EM_lumpedMat, q, s, isub;
41 type_t funcspace;
42 index_t color,*row_index=NULL;
43 __const double *D_p=NULL;
44 double *S=NULL, *EM_lumpedMat=NULL, *Vol=NULL, *lumpedMat_p=NULL;
45 register double rtmp;
46 size_t len_EM_lumpedMat_size;
47 register double m_t=0., diagS=0., rtmp2=0.;
48
49 Finley_resetError();
50
51 if (nodes==NULL || elements==NULL) return;
52 if (isEmpty(lumpedMat) || isEmpty(D)) return;
53 if (isEmpty(lumpedMat) && !isEmpty(D)) {
54 Finley_setError(TYPE_ERROR,"Finley_Assemble_LumpedSystem: coefficients are non-zero but no lumped matrix is given.");
55 return;
56 }
57 funcspace=getFunctionSpaceType(D);
58 /* check if all function spaces are the same */
59 if (funcspace==FINLEY_ELEMENTS) {
60 reducedIntegrationOrder=FALSE;
61 } else if (funcspace==FINLEY_FACE_ELEMENTS) {
62 reducedIntegrationOrder=FALSE;
63 } else if (funcspace==FINLEY_REDUCED_ELEMENTS) {
64 reducedIntegrationOrder=TRUE;
65 } else if (funcspace==FINLEY_REDUCED_FACE_ELEMENTS) {
66 reducedIntegrationOrder=TRUE;
67 } else {
68 Finley_setError(TYPE_ERROR,"Finley_Assemble_LumpedSystem: assemblage failed because of illegal function space.");
69 }
70 if (! Finley_noError()) return;
71
72 /* set all parameters in p*/
73 Finley_Assemble_getAssembleParameters(nodes,elements,NULL,lumpedMat, reducedIntegrationOrder, &p);
74 if (! Finley_noError()) return;
75
76 /* check if all function spaces are the same */
77 if (! numSamplesEqual(D,p.numQuadTotal,elements->numElements) ) {
78 sprintf(error_msg,"Finley_Assemble_LumpedSystem: sample points of coefficient D don't match (%d,%d)",p.numQuadSub,elements->numElements);
79 Finley_setError(TYPE_ERROR,error_msg);
80 }
81
82 /* check the dimensions: */
83 if (p.numEqu==1) {
84 if (!isEmpty(D)) {
85 if (!isDataPointShapeEqual(D,0,dimensions)) {
86 Finley_setError(TYPE_ERROR,"Finley_Assemble_LumpedSystem: coefficient D, rank 0 expected.");
87 }
88
89 }
90 } else {
91 if (!isEmpty(D)) {
92 dimensions[0]=p.numEqu;
93 if (!isDataPointShapeEqual(D,1,dimensions)) {
94 sprintf(error_msg,"Finley_Assemble_LumpedSystem: coefficient D, expected shape (%d,)",dimensions[0]);
95 Finley_setError(TYPE_ERROR,error_msg);
96 }
97 }
98 }
99 if (Finley_noError()) {
100 requireWrite(lumpedMat);
101 lumpedMat_p=getSampleDataRW(lumpedMat,0);
102 len_EM_lumpedMat=p.row_numShapesTotal*p.numEqu;
103 len_EM_lumpedMat_size=len_EM_lumpedMat*sizeof(double);
104
105 expandedD=isExpanded(D);
106 S=p.row_jac->BasisFunctions->S;
107
108 #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp, diagS, m_t, isub, rtmp2)
109 {
110 EM_lumpedMat=THREAD_MEMALLOC(len_EM_lumpedMat,double);
111 row_index=THREAD_MEMALLOC(p.row_numShapesTotal,index_t);
112 if ( !Finley_checkPtr(EM_lumpedMat) && !Finley_checkPtr(row_index) ) {
113 if (p.numEqu == 1) { /* single equation */
114 if (expandedD) { /* with expanded D */
115 for (color=elements->minColor;color<=elements->maxColor;color++) {
116 /* open loop over all elements: */
117 #pragma omp for private(e) schedule(static)
118 for(e=0;e<elements->numElements;e++){
119 if (elements->Color[e]==color) {
120 for (isub=0; isub<p.numSub; isub++) {
121 Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);
122 D_p=getSampleDataRO(D,e);
123 if (useHRZ) {
124
125 m_t=0; /* mass of the element: m_t */
126 #pragma ivdep
127 for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q]*D_p[INDEX2(q, isub,p.numQuadSub) ];
128
129 diagS=0; /* diagonal sum: S */
130 for (s=0;s<p.row_numShapes;s++) {
131 rtmp=0;
132 #pragma ivdep
133 for (q=0;q<p.numQuadSub;q++) {
134 rtmp2=S[INDEX2(s,q,p.row_numShapes)];
135 rtmp+=Vol[q]*D_p[INDEX2(q, isub,p.numQuadSub)]*rtmp2*rtmp2;
136 }
137 EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;
138 diagS+=rtmp;
139 }
140 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
141 rtmp=m_t/diagS;
142 #pragma ivdep
143 for (s=0;s<p.row_numShapes;s++) EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=rtmp;
144
145 } else { /* row-sum lumping */
146 for (s=0;s<p.row_numShapes;s++) {
147 rtmp=0;
148 #pragma ivdep
149 for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*D_p[INDEX2(q, isub,p.numQuadSub)];
150 EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;
151 }
152
153 }
154 for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[INDEX2(q,isub,p.row_numShapesTotal)],e,p.NN)]];
155 Finley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
156 } /* end of isub loop */
157 } /* end color check */
158 } /* end element loop */
159 } /* end color loop */
160
161
162 } else { /* with constant D */
163
164 for (color=elements->minColor;color<=elements->maxColor;color++) {
165 /* open loop over all elements: */
166 #pragma omp for private(e) schedule(static)
167 for(e=0;e<elements->numElements;e++){
168 if (elements->Color[e]==color) {
169 for (isub=0; isub<p.numSub; isub++) {
170 Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);
171 D_p=getSampleDataRO(D,e);
172 if (useHRZ) { /* HRZ lumping */
173 m_t=0; /* mass of the element: m_t */
174 #pragma ivdep
175 for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q];
176 diagS=0; /* diagonal sum: S */
177 for (s=0;s<p.row_numShapes;s++) {
178 rtmp=0;
179 #pragma ivdep
180 for (q=0;q<p.numQuadSub;q++){
181 rtmp2=S[INDEX2(s,q,p.row_numShapes)];
182 rtmp+=Vol[q]*rtmp2*rtmp2;
183 }
184 EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;
185 diagS+=rtmp;
186 }
187 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
188 rtmp=m_t/diagS*D_p[0];
189 #pragma ivdep
190 for (s=0;s<p.row_numShapes;s++) EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=rtmp;
191 } else { /* row-sum lumping */
192 for (s=0;s<p.row_numShapes;s++) {
193 rtmp=0;
194 #pragma ivdep
195 for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];
196 EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp*D_p[0];
197 }
198 }
199 for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[INDEX2(q,isub,p.row_numShapesTotal)],e,p.NN)]];
200 Finley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
201 } /* end of isub loop */
202 } /* end color check */
203 } /* end element loop */
204 } /* end color loop */
205
206 }
207 } else { /* system of equation */
208 if (expandedD) { /* with expanded D */
209 for (color=elements->minColor;color<=elements->maxColor;color++) {
210 /* open loop over all elements: */
211 #pragma omp for private(e) schedule(static)
212 for(e=0;e<elements->numElements;e++){
213 if (elements->Color[e]==color) {
214 for (isub=0; isub<p.numSub; isub++) {
215 Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);
216 D_p=getSampleDataRO(D,e);
217
218 if (useHRZ) { /* HRZ lumping */
219 for (k=0;k<p.numEqu;k++) {
220 m_t=0; /* mass of the element: m_t */
221 #pragma ivdep
222 for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q]*D_p[INDEX3(k,q,isub,p.numEqu,p.numQuadSub)];
223
224 diagS=0; /* diagonal sum: S */
225 for (s=0;s<p.row_numShapes;s++) {
226 rtmp=0;
227 #pragma ivdep
228 for (q=0;q<p.numQuadSub;q++) {
229 rtmp2=S[INDEX2(s,q,p.row_numShapes)];
230 rtmp+=Vol[q]*D_p[INDEX3(k,q,isub,p.numEqu,p.numQuadSub)]*rtmp2*rtmp2;
231 }
232 EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
233 diagS+=rtmp;
234 }
235 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
236 rtmp=m_t/diagS;
237 #pragma ivdep
238 for (s=0;s<p.row_numShapes;s++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=rtmp;
239 }
240 } else { /* row-sum lumping */
241 for (s=0;s<p.row_numShapes;s++) {
242 for (k=0;k<p.numEqu;k++) {
243 rtmp=0.;
244 #pragma ivdep
245 for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*D_p[INDEX3(k,q,isub,p.numEqu,p.numQuadSub)];
246 EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
247 }
248 }
249 }
250 for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[INDEX2(q,isub,p.row_numShapesTotal)],e,p.NN)]];
251 Finley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
252 } /* end of isub loop */
253 } /* end color check */
254 } /* end element loop */
255 } /* end color loop */
256 } else { /* with constant D */
257 for (color=elements->minColor;color<=elements->maxColor;color++) {
258 /* open loop over all elements: */
259 #pragma omp for private(e) schedule(static)
260 for(e=0;e<elements->numElements;e++){
261 if (elements->Color[e]==color) {
262 for (isub=0; isub<p.numSub; isub++) {
263 Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);
264 D_p=getSampleDataRO(D,e);
265
266 if (useHRZ) { /* HRZ lumping */
267 m_t=0; /* mass of the element: m_t */
268 #pragma ivdep
269 for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q];
270 diagS=0; /* diagonal sum: S */
271 for (s=0;s<p.row_numShapes;s++) {
272 rtmp=0;
273 #pragma ivdep
274 for (q=0;q<p.numQuadSub;q++) {
275 rtmp2=S[INDEX2(s,q,p.row_numShapes)];
276 rtmp+=Vol[q]*rtmp2*rtmp2;
277 }
278 #pragma ivdep
279 for (k=0;k<p.numEqu;k++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
280 diagS+=rtmp;
281 }
282
283 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
284 rtmp=m_t/diagS;
285 for (s=0;s<p.row_numShapes;s++) {
286 #pragma ivdep
287 for (k=0;k<p.numEqu;k++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=rtmp*D_p[k];
288 }
289 } else { /* row-sum lumping */
290 for (s=0;s<p.row_numShapes;s++) {
291 for (k=0;k<p.numEqu;k++) {
292 rtmp=0.;
293 #pragma ivdep
294 for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];
295 EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp*D_p[k];
296 }
297 }
298 }
299 for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[INDEX2(q,isub,p.row_numShapesTotal)],e,p.NN)]];
300 Finley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
301 } /* end of isub loop */
302 } /* end color check */
303 } /* end element loop */
304 } /* end color loop */
305 }
306 }
307 } /* end of pointer check */
308 THREAD_MEMFREE(EM_lumpedMat);
309 THREAD_MEMFREE(row_index);
310 } /* end parallel region */
311 }
312 }

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26