/[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 2770 - (hide annotations)
Wed Nov 25 01:24:51 2009 UTC (9 years, 3 months ago) by jfenwick
File MIME type: text/plain
File size: 15182 byte(s)
Removed buffer implementation of Lazy.
Removed supporting Alloc/Free Sample buffer calls.
1 jgs 82
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4 jfenwick 2548 * Copyright (c) 2003-2009 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 ksteube 1221 /* Disabled until the tests pass */
33     /* #define NEW_LUMPING */ /* */
34 btully 1220
35 jgs 82 /**************************************************************/
36    
37 gross 1204 void Finley_Assemble_LumpedSystem(Finley_NodeFile* nodes,Finley_ElementFile* elements, escriptDataC* lumpedMat, escriptDataC* D)
38     {
39 jgs 82
40 gross 1204 bool_t reducedIntegrationOrder=FALSE, expandedD;
41 jgs 150 char error_msg[LenErrorMsg_MAX];
42 gross 798 Assemble_Parameters p;
43 gross 2748 dim_t dimensions[ESCRIPT_MAX_DATA_RANK], k, e, len_EM_lumpedMat, q, s, isub;
44 gross 1028 type_t funcspace;
45 gross 1204 index_t color,*row_index=NULL;
46 jfenwick 2271 __const double *D_p=NULL;
47     double *S=NULL, *EM_lumpedMat=NULL, *Vol=NULL, *lumpedMat_p=NULL;
48 phornby 1628 register double rtmp;
49 gross 1204 size_t len_EM_lumpedMat_size;
50 caltinay 2750 #ifdef NEW_LUMPING
51 gross 2748 register double m_t=0., diagS=0.;
52 caltinay 2750 #endif
53 gross 2748
54 jgs 150 Finley_resetError();
55 jgs 82
56     if (nodes==NULL || elements==NULL) return;
57 gross 1204 if (isEmpty(lumpedMat) || isEmpty(D)) return;
58     if (isEmpty(lumpedMat) && !isEmpty(D)) {
59     Finley_setError(TYPE_ERROR,"Assemble_LumpedSystem: coefficients are non-zero but no lumped matrix is given.");
60     return;
61 gross 798 }
62 gross 1204 funcspace=getFunctionSpaceType(D);
63 jgs 82 /* check if all function spaces are the same */
64 gross 798 if (funcspace==FINLEY_ELEMENTS) {
65     reducedIntegrationOrder=FALSE;
66     } else if (funcspace==FINLEY_FACE_ELEMENTS) {
67     reducedIntegrationOrder=FALSE;
68 gross 1072 } else if (funcspace==FINLEY_REDUCED_ELEMENTS) {
69     reducedIntegrationOrder=TRUE;
70     } else if (funcspace==FINLEY_REDUCED_FACE_ELEMENTS) {
71     reducedIntegrationOrder=TRUE;
72 gross 798 } else {
73 gross 1204 Finley_setError(TYPE_ERROR,"Assemble_LumpedSystem: assemblage failed because of illegal function space.");
74 gross 798 }
75     if (! Finley_noError()) return;
76 jgs 82
77 gross 798 /* set all parameters in p*/
78 gross 1204 Assemble_getAssembleParameters(nodes,elements,NULL,lumpedMat, reducedIntegrationOrder, &p);
79 gross 798 if (! Finley_noError()) return;
80    
81     /* check if all function spaces are the same */
82    
83 gross 2748 if (! numSamplesEqual(D,p.numQuadSub,elements->numElements) ) {
84     sprintf(error_msg,"Assemble_LumpedSystem: sample points of coefficient D don't match (%d,%d)",p.numQuadSub,elements->numElements);
85 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
86 jgs 82 }
87    
88     /* check the dimensions: */
89    
90 gross 1204 if (p.numEqu==1) {
91 jgs 82 if (!isEmpty(D)) {
92     if (!isDataPointShapeEqual(D,0,dimensions)) {
93 gross 1204 Finley_setError(TYPE_ERROR,"Assemble_LumpedSystem: coefficient D, rank 0 expected.");
94 jgs 82 }
95 gross 1204
96 jgs 82 }
97     } else {
98     if (!isEmpty(D)) {
99     dimensions[0]=p.numEqu;
100 gross 1204 if (!isDataPointShapeEqual(D,1,dimensions)) {
101     sprintf(error_msg,"Assemble_LumpedSystem: coefficient D, expected shape (%d,)",dimensions[0]);
102 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
103 jgs 82 }
104     }
105     }
106 gross 1204
107 jgs 150 if (Finley_noError()) {
108 jfenwick 2271 requireWrite(lumpedMat);
109     lumpedMat_p=getSampleDataRW(lumpedMat,0);
110 gross 2748 len_EM_lumpedMat=p.row_numShapesTotal*p.numEqu;
111 gross 1204 len_EM_lumpedMat_size=len_EM_lumpedMat*sizeof(double);
112 gross 2748
113 gross 1204 expandedD=isExpanded(D);
114 gross 2748 S=p.row_jac->BasisFunctions->S;
115 caltinay 2750
116     #ifdef NEW_LUMPING
117 gross 2748 #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp, diagS, m_t)
118 caltinay 2750 #else
119     #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp)
120     #endif
121 gross 1204 {
122     EM_lumpedMat=THREAD_MEMALLOC(len_EM_lumpedMat,double);
123 gross 2748 row_index=THREAD_MEMALLOC(p.row_numShapesTotal,index_t);
124 gross 1204 if ( !Finley_checkPtr(EM_lumpedMat) && !Finley_checkPtr(row_index) ) {
125     if (p.numEqu == 1) {
126     if (expandedD) {
127 gross 2748
128 gross 1204 for (color=elements->minColor;color<=elements->maxColor;color++) {
129     /* open loop over all elements: */
130     #pragma omp for private(e) schedule(static)
131     for(e=0;e<elements->numElements;e++){
132 gross 2748
133     if (elements->Color[e]==color) {
134     for (isub=0; isub<p.numSub; isub++) {
135     Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);
136     memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
137    
138 jfenwick 2770 D_p=getSampleDataRO(D,e);
139 gross 2748 #ifdef NEW_LUMPING /* HRZ lumping */
140     m_t=0; /* mass of the element: m_t */
141     for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q]*D_p[INDEX2(q, isub,p.numQuadSub) ];
142    
143     diagS=0; /* diagonal sum: S */
144     for (s=0;s<p.row_numShapes;s++) {
145     rtmp=0;
146     for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*D_p[INDEX2(q, isub,p.numQuadSub)]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];
147     EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;
148     diagS+=rtmp;
149     }
150     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
151     rtmp=m_t/diagS;
152     for (s=0;s<p.row_numShapes;s++) EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=rtmp;
153    
154     #else /* row-sum lumping */
155     for (s=0;s<p.row_numShapes;s++) {
156     rtmp=0;
157     for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*D_p[INDEX2(q, isub,p.numQuadSub)];
158     EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
159     }
160     #endif
161     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)]];
162     Finley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
163     } /* end of isub loop */
164 gross 1204 } /* end color check */
165 gross 2748
166 gross 1204 } /* end element loop */
167     } /* end color loop */
168 gross 2748 } else {
169 gross 1204 for (color=elements->minColor;color<=elements->maxColor;color++) {
170 gross 2748 /* open loop over all elements: */
171     #pragma omp for private(e) schedule(static)
172     for(e=0;e<elements->numElements;e++){
173    
174     if (elements->Color[e]==color) {
175     for (isub=0; isub<p.numSub; isub++) {
176    
177     Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);
178     memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
179    
180 jfenwick 2770 D_p=getSampleDataRO(D,e);
181 gross 2748 #ifdef NEW_LUMPING /* HRZ lumping */
182     m_t=0; /* mass of the element: m_t */
183     for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q];
184     m_t*=D_p[0];
185    
186     diagS=0; /* diagonal sum: S */
187     for (s=0;s<p.row_numShapes;s++) {
188     rtmp=0;
189     for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];
190     rtmp*=D_p[0];
191     EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;
192     diagS+=rtmp;
193     }
194     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
195     rtmp=m_t/diagS;
196     for (s=0;s<p.row_numShapes;s++) EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=rtmp;
197    
198     #else /* row-sum lumping */
199     for (s=0;s<p.row_numShapes;s++) {
200     rtmp=0;
201     for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];
202     EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp*D_p[0];
203     }
204     #endif
205     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)]];
206     Finley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
207     } /* end of isub loop */
208     } /* end color check */
209 gross 1204 } /* end element loop */
210 gross 2748 } /* end color loop */
211    
212 gross 1204 }
213 gross 798 } else {
214 gross 1204 if (expandedD) {
215     for (color=elements->minColor;color<=elements->maxColor;color++) {
216     /* open loop over all elements: */
217     #pragma omp for private(e) schedule(static)
218     for(e=0;e<elements->numElements;e++){
219     if (elements->Color[e]==color) {
220 gross 2748 for (isub=0; isub<p.numSub; isub++) {
221     Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);
222     memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
223 jfenwick 2770 D_p=getSampleDataRO(D,e);
224 gross 2748
225     #ifdef NEW_LUMPING /* HRZ lumping */
226     for (k=0;k<p.numEqu;k++) {
227     m_t=0; /* mass of the element: m_t */
228     for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q]*D_p[INDEX3(k,q,isub,p.numEqu,p.numQuadSub)];
229    
230     diagS=0; /* diagonal sum: S */
231     for (s=0;s<p.row_numShapes;s++) {
232     rtmp=0;
233     for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*D_p[INDEX3(k,q,isub,p.numEqu,p.numQuadSub)]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];
234     EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
235     diagS+=rtmp;
236     }
237     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
238     rtmp=m_t/diagS;
239     for (s=0;s<p.row_numShapes;s++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=rtmp;
240     }
241     #else /* row-sum lumping */
242     for (s=0;s<p.row_numShapes;s++) {
243     for (k=0;k<p.numEqu;k++) {
244     rtmp=0.;
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 btully 1220 }
248 gross 2748 }
249     #endif
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 gross 1204 } /* end color check */
254     } /* end element loop */
255     } /* end color loop */
256     } else {
257     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 2748 for (isub=0; isub<p.numSub; isub++) {
263     Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);
264     memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
265 jfenwick 2770 D_p=getSampleDataRO(D,e);
266 gross 2748
267     #ifdef NEW_LUMPING /* HRZ lumping */
268     for (k=0;k<p.numEqu;k++) {
269     m_t=0; /* mass of the element: m_t */
270     for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q]*D_p[INDEX3(k,q,isub,p.numEqu,p.numQuadSub)];
271     m_t*=D_p[k];
272     diagS=0; /* diagonal sum: S */
273     for (s=0;s<p.row_numShapes;s++) {
274     rtmp=0;
275     for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];
276     rtmp*=D_p[k];
277     EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
278     diagS+=rtmp;
279     }
280     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
281     rtmp=m_t/diagS;
282     for (s=0;s<p.row_numShapes;s++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=rtmp;
283     }
284     #else /* row-sum lumping */
285     for (s=0;s<p.row_numShapes;s++) {
286     for (k=0;k<p.numEqu;k++) {
287     rtmp=0.;
288     for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];
289     EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp*D_p[k];
290     }
291     }
292     #endif
293     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)]];
294     Finley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
295     } /* end of isub loop */
296 gross 1204 } /* end color check */
297     } /* end element loop */
298     } /* end color loop */
299     }
300 gross 798 }
301 gross 1204 } /* end of pointer check */
302     THREAD_MEMFREE(EM_lumpedMat);
303     THREAD_MEMFREE(row_index);
304     } /* end parallel region */
305 jgs 82 }
306     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26