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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26