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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3379 - (hide annotations)
Wed Nov 24 04:48:49 2010 UTC (8 years, 4 months ago) by gross
File MIME type: text/plain
File size: 15065 byte(s)
some clarification on lumping
1 jgs 82
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4 jfenwick 2881 * Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1811 * 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 ksteube 1312
14 ksteube 1811
15 jgs 82 /**************************************************************/
16    
17 gross 1204 /* assembles the mass matrix in lumped form */
18 jgs 82
19 gross 1204 /* The coefficient D has to be defined on the integration points or not present. */
20 jgs 82
21 gross 1204 /* lumpedMat has to be initialized before the routine is called. */
22 jgs 82
23     /**************************************************************/
24    
25     #include "Assemble.h"
26     #include "Util.h"
27     #ifdef _OPENMP
28     #include <omp.h>
29     #endif
30    
31    
32     /**************************************************************/
33    
34 gross 3379 void Finley_Assemble_LumpedSystem(Finley_NodeFile* nodes,Finley_ElementFile* elements, escriptDataC* lumpedMat, escriptDataC* D, const bool_t useHRZ)
35 gross 1204 {
36 jgs 82
37 gross 1204 bool_t reducedIntegrationOrder=FALSE, expandedD;
38 jgs 150 char error_msg[LenErrorMsg_MAX];
39 jfenwick 3259 Finley_Assemble_Parameters p;
40 gross 2748 dim_t dimensions[ESCRIPT_MAX_DATA_RANK], k, e, len_EM_lumpedMat, q, s, isub;
41 gross 1028 type_t funcspace;
42 gross 1204 index_t color,*row_index=NULL;
43 jfenwick 2271 __const double *D_p=NULL;
44     double *S=NULL, *EM_lumpedMat=NULL, *Vol=NULL, *lumpedMat_p=NULL;
45 phornby 1628 register double rtmp;
46 gross 1204 size_t len_EM_lumpedMat_size;
47 gross 3379 register double m_t=0., diagS=0., rtmp2=0.;
48 gross 2989
49 jgs 150 Finley_resetError();
50 jgs 82
51     if (nodes==NULL || elements==NULL) return;
52 gross 1204 if (isEmpty(lumpedMat) || isEmpty(D)) return;
53     if (isEmpty(lumpedMat) && !isEmpty(D)) {
54 jfenwick 3259 Finley_setError(TYPE_ERROR,"Finley_Assemble_LumpedSystem: coefficients are non-zero but no lumped matrix is given.");
55 gross 1204 return;
56 gross 798 }
57 gross 1204 funcspace=getFunctionSpaceType(D);
58 jgs 82 /* check if all function spaces are the same */
59 gross 798 if (funcspace==FINLEY_ELEMENTS) {
60     reducedIntegrationOrder=FALSE;
61     } else if (funcspace==FINLEY_FACE_ELEMENTS) {
62     reducedIntegrationOrder=FALSE;
63 gross 1072 } else if (funcspace==FINLEY_REDUCED_ELEMENTS) {
64     reducedIntegrationOrder=TRUE;
65     } else if (funcspace==FINLEY_REDUCED_FACE_ELEMENTS) {
66     reducedIntegrationOrder=TRUE;
67 gross 798 } else {
68 jfenwick 3259 Finley_setError(TYPE_ERROR,"Finley_Assemble_LumpedSystem: assemblage failed because of illegal function space.");
69 gross 798 }
70     if (! Finley_noError()) return;
71 jgs 82
72 gross 798 /* set all parameters in p*/
73 jfenwick 3259 Finley_Assemble_getAssembleParameters(nodes,elements,NULL,lumpedMat, reducedIntegrationOrder, &p);
74 gross 798 if (! Finley_noError()) return;
75 gross 2989
76 gross 798 /* check if all function spaces are the same */
77 gross 2989 if (! numSamplesEqual(D,p.numQuadTotal,elements->numElements) ) {
78 jfenwick 3259 sprintf(error_msg,"Finley_Assemble_LumpedSystem: sample points of coefficient D don't match (%d,%d)",p.numQuadSub,elements->numElements);
79 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
80 jgs 82 }
81    
82     /* check the dimensions: */
83 gross 1204 if (p.numEqu==1) {
84 jgs 82 if (!isEmpty(D)) {
85     if (!isDataPointShapeEqual(D,0,dimensions)) {
86 jfenwick 3259 Finley_setError(TYPE_ERROR,"Finley_Assemble_LumpedSystem: coefficient D, rank 0 expected.");
87 jgs 82 }
88 gross 1204
89 jgs 82 }
90     } else {
91     if (!isEmpty(D)) {
92     dimensions[0]=p.numEqu;
93 gross 1204 if (!isDataPointShapeEqual(D,1,dimensions)) {
94 jfenwick 3259 sprintf(error_msg,"Finley_Assemble_LumpedSystem: coefficient D, expected shape (%d,)",dimensions[0]);
95 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
96 jgs 82 }
97     }
98     }
99 jgs 150 if (Finley_noError()) {
100 jfenwick 2271 requireWrite(lumpedMat);
101     lumpedMat_p=getSampleDataRW(lumpedMat,0);
102 gross 2748 len_EM_lumpedMat=p.row_numShapesTotal*p.numEqu;
103 gross 1204 len_EM_lumpedMat_size=len_EM_lumpedMat*sizeof(double);
104 gross 2748
105 gross 1204 expandedD=isExpanded(D);
106 gross 2748 S=p.row_jac->BasisFunctions->S;
107 gross 2989
108 gross 3379 #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp, diagS, m_t, isub, rtmp2)
109 gross 1204 {
110     EM_lumpedMat=THREAD_MEMALLOC(len_EM_lumpedMat,double);
111 gross 2748 row_index=THREAD_MEMALLOC(p.row_numShapesTotal,index_t);
112 gross 1204 if ( !Finley_checkPtr(EM_lumpedMat) && !Finley_checkPtr(row_index) ) {
113 gross 2988 if (p.numEqu == 1) { /* single equation */
114     if (expandedD) { /* with expanded D */
115 gross 1204 for (color=elements->minColor;color<=elements->maxColor;color++) {
116     /* open loop over all elements: */
117     #pragma omp for private(e) schedule(static)
118 gross 2988 for(e=0;e<elements->numElements;e++){
119     if (elements->Color[e]==color) {
120 gross 2748 for (isub=0; isub<p.numSub; isub++) {
121 gross 2989 Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);
122 jfenwick 2770 D_p=getSampleDataRO(D,e);
123 gross 3379 if (useHRZ) {
124 gross 2988
125 gross 2748 m_t=0; /* mass of the element: m_t */
126 gross 3379 #pragma ivdep
127 gross 2748 for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q]*D_p[INDEX2(q, isub,p.numQuadSub) ];
128 gross 2988
129 gross 2748 diagS=0; /* diagonal sum: S */
130     for (s=0;s<p.row_numShapes;s++) {
131     rtmp=0;
132 gross 3379 #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 gross 2748 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 gross 2988 rtmp=m_t/diagS;
142 gross 3379 #pragma ivdep
143 gross 2748 for (s=0;s<p.row_numShapes;s++) EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=rtmp;
144    
145 gross 3379 } else { /* row-sum lumping */
146 gross 2748 for (s=0;s<p.row_numShapes;s++) {
147     rtmp=0;
148 gross 3379 #pragma ivdep
149 gross 2748 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 gross 2989 EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;
151 gross 2748 }
152 gross 2989
153 gross 3379 }
154 gross 2748 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 gross 3379 } /* end color check */
158 gross 1204 } /* end element loop */
159     } /* end color loop */
160 gross 3379
161    
162 gross 2988 } else { /* with constant D */
163 gross 2990
164 gross 1204 for (color=elements->minColor;color<=elements->maxColor;color++) {
165 gross 2748 /* open loop over all elements: */
166     #pragma omp for private(e) schedule(static)
167     for(e=0;e<elements->numElements;e++){
168 gross 2988 if (elements->Color[e]==color) {
169 gross 2748 for (isub=0; isub<p.numSub; isub++) {
170 gross 2988 Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);
171 jfenwick 2770 D_p=getSampleDataRO(D,e);
172 gross 3379 if (useHRZ) { /* HRZ lumping */
173 gross 2748 m_t=0; /* mass of the element: m_t */
174 gross 3379 #pragma ivdep
175 gross 2748 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 gross 3379 #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 gross 2748 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 gross 2988 rtmp=m_t/diagS*D_p[0];
189 gross 3379 #pragma ivdep
190 gross 2748 for (s=0;s<p.row_numShapes;s++) EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=rtmp;
191 gross 3379 } else { /* row-sum lumping */
192 gross 2748 for (s=0;s<p.row_numShapes;s++) {
193     rtmp=0;
194 gross 3379 #pragma ivdep
195 gross 2748 for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];
196 gross 2988 EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp*D_p[0];
197 gross 2748 }
198 gross 3379 }
199 gross 2748 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 gross 1204 } /* end element loop */
204 gross 2748 } /* end color loop */
205    
206 gross 1204 }
207 gross 2988 } else { /* system of equation */
208     if (expandedD) { /* with expanded D */
209 gross 1204 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 gross 2988 for (isub=0; isub<p.numSub; isub++) {
215 gross 2748 Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);
216 gross 2988 D_p=getSampleDataRO(D,e);
217    
218 gross 3379 if (useHRZ) { /* HRZ lumping */
219 gross 2748 for (k=0;k<p.numEqu;k++) {
220     m_t=0; /* mass of the element: m_t */
221 gross 3379 #pragma ivdep
222 gross 2748 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 gross 3379 #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 gross 2748 EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
233     diagS+=rtmp;
234 gross 2988 }
235     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
236     rtmp=m_t/diagS;
237 gross 3379 #pragma ivdep
238 gross 2988 for (s=0;s<p.row_numShapes;s++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=rtmp;
239     }
240 gross 3379 } else { /* row-sum lumping */
241 gross 2988 for (s=0;s<p.row_numShapes;s++) {
242     for (k=0;k<p.numEqu;k++) {
243     rtmp=0.;
244 gross 3379 #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 gross 2988 }
248     }
249 gross 3379 }
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 gross 2988 } /* end of isub loop */
253 gross 1204 } /* end color check */
254     } /* end element loop */
255     } /* end color loop */
256 gross 2988 } else { /* with constant D */
257 gross 1204 for (color=elements->minColor;color<=elements->maxColor;color++) {
258 gross 2748 /* open loop over all elements: */
259 gross 1204 #pragma omp for private(e) schedule(static)
260     for(e=0;e<elements->numElements;e++){
261     if (elements->Color[e]==color) {
262 gross 2988 for (isub=0; isub<p.numSub; isub++) {
263 gross 2748 Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);
264 jfenwick 2770 D_p=getSampleDataRO(D,e);
265 gross 2748
266 gross 3379 if (useHRZ) { /* HRZ lumping */
267 gross 2748 m_t=0; /* mass of the element: m_t */
268 gross 3379 #pragma ivdep
269 gross 2988 for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q];
270 gross 2748 diagS=0; /* diagonal sum: S */
271     for (s=0;s<p.row_numShapes;s++) {
272     rtmp=0;
273 gross 3379 #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 gross 2988 for (k=0;k<p.numEqu;k++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
280 gross 2748 diagS+=rtmp;
281     }
282 gross 3379
283 gross 2988 /* 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 gross 3379 #pragma ivdep
287 gross 2988 for (k=0;k<p.numEqu;k++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=rtmp*D_p[k];
288     }
289 gross 3379 } else { /* row-sum lumping */
290 gross 2748 for (s=0;s<p.row_numShapes;s++) {
291 gross 2988 for (k=0;k<p.numEqu;k++) {
292 gross 2748 rtmp=0.;
293 gross 3379 #pragma ivdep
294 gross 2748 for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];
295 gross 2989 EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp*D_p[k];
296 gross 2748 }
297 gross 2988 }
298 gross 3379 }
299 gross 2988 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 gross 1204 } /* end color check */
303     } /* end element loop */
304     } /* end color loop */
305     }
306 gross 798 }
307 gross 1204 } /* end of pointer check */
308     THREAD_MEMFREE(EM_lumpedMat);
309     THREAD_MEMFREE(row_index);
310     } /* end parallel region */
311 jgs 82 }
312     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26