/[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 1221 - (hide annotations)
Fri Aug 3 00:27:20 2007 UTC (11 years, 8 months ago) by ksteube
File MIME type: text/plain
File size: 15940 byte(s)
Disabled HRZ lumping until the tests pass

1 jgs 150 /*
2 elspeth 616 ************************************************************
3 gross 1204 * Copyright 2006,2007 by ACcENULLNULL MNRF *
4 elspeth 616 * *
5     * http://www.access.edu.au *
6     * Primary Business: Queensland, Australia *
7 gross 1204 * Licensed under the Open NULLoftware License version 3.0 *
8 elspeth 616 * http://www.opensource.org/licenses/osl-3.0.php *
9     * *
10     ************************************************************
11 jgs 150 */
12 jgs 82
13 jgs 150
14 jgs 82 /**************************************************************/
15    
16 gross 1204 /* assembles the mass matrix in lumped form */
17 jgs 82
18 gross 1204 /* The coefficient D has to be defined on the integration points or not present. */
19 jgs 82
20 gross 1204 /* lumpedMat has to be initialized before the routine is called. */
21 jgs 82
22     /**************************************************************/
23    
24 jgs 150 /* Author: gross@access.edu.au */
25 btully 1220 /* Version: $Id$ */
26 jgs 82
27     /**************************************************************/
28    
29     #include "Assemble.h"
30     #include "Util.h"
31     #ifdef _OPENMP
32     #include <omp.h>
33     #endif
34    
35    
36 ksteube 1221 /* Disabled until the tests pass */
37     /* #define NEW_LUMPING */ /* */
38 btully 1220
39 jgs 82 /**************************************************************/
40    
41 gross 1204 void Finley_Assemble_LumpedSystem(Finley_NodeFile* nodes,Finley_ElementFile* elements, escriptDataC* lumpedMat, escriptDataC* D)
42     {
43 jgs 82
44 gross 1204 bool_t reducedIntegrationOrder=FALSE, expandedD;
45 jgs 150 char error_msg[LenErrorMsg_MAX];
46 gross 798 Assemble_Parameters p;
47 jgs 82 double time0;
48 gross 1204 dim_t dimensions[ESCRIPT_MAX_DATA_RANK], k, e, len_EM_lumpedMat, q, s;
49 gross 1028 type_t funcspace;
50 gross 1204 index_t color,*row_index=NULL;
51     double *S=NULL, *EM_lumpedMat=NULL, *Vol=NULL, *D_p=NULL, *lumpedMat_p=NULL;
52 btully 1220 register double rtmp, m_t, diagS;
53 gross 1204 size_t len_EM_lumpedMat_size;
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 btully 1220 #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp, m_t, diagS)
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     #ifndef PASO_MPI
123     for (color=elements->minColor;color<=elements->maxColor;color++) {
124     /* open loop over all elements: */
125     #pragma omp for private(e) schedule(static)
126     for(e=0;e<elements->numElements;e++){
127     if (elements->Color[e]==color) {
128     #else
129     {
130     for(e=0;e<elements->numElements;e++) {
131     {
132     #endif
133     Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
134     memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
135     D_p=getSampleData(D,e);
136 btully 1220 #ifdef NEW_LUMPING /* HRZ lumping */
137     /*
138     * Number of PDEs: 1
139     * D_p varies over element: True
140     */
141     m_t=0; /* mass of the element: m_t */
142     for (q=0;q<p.numQuad;q++) {
143     m_t+=Vol[q]*D_p[q];
144     }
145     diagS=0; /* diagonal sum: S */
146 gross 1204 for (s=0;s<p.row_NS;s++) {
147     rtmp=0;
148 btully 1220 for (q=0;q<p.numQuad;q++) {
149     rtmp+=Vol[q]*D_p[q]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
150     }
151     EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
152     diagS+=EM_lumpedMat[INDEX2(0,s,p.numEqu)];
153     }
154     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
155     for (s=0;s<p.row_NS;s++) {
156     EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=m_t/diagS;
157     }
158     #else /* row-sum lumping */
159     for (s=0;s<p.row_NS;s++) {
160     rtmp=0;
161 gross 1204 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[q];
162     EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
163     }
164 btully 1220 #endif
165 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)]];
166     Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
167     } /* end color check */
168     } /* end element loop */
169     } /* end color loop */
170     } else {
171     #ifndef PASO_MPI
172     for (color=elements->minColor;color<=elements->maxColor;color++) {
173     /* open loop over all elements: */
174     #pragma omp for private(e) schedule(static)
175     for(e=0;e<elements->numElements;e++){
176     if (elements->Color[e]==color) {
177     #else
178     {
179     for(e=0;e<elements->numElements;e++) {
180     {
181     #endif
182 btully 1220 Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
183     memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
184     D_p=getSampleData(D,e);
185     #ifdef NEW_LUMPING /* HRZ lumping */
186     /*
187     * Number of PDEs: 1
188     * D_p varies over element: False
189     */
190     m_t=0; /* mass of the element: m_t */
191     for (q=0;q<p.numQuad;q++) {
192     m_t+=Vol[q]*D_p[0];
193     }
194     diagS=0; /* diagonal sum: S */
195     for (s=0;s<p.row_NS;s++) {
196     rtmp=0;
197     for (q=0;q<p.numQuad;q++) {
198     rtmp+=Vol[q]*D_p[0]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
199     }
200     EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
201     diagS+=EM_lumpedMat[INDEX2(0,s,p.numEqu)];
202     }
203     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
204     for (s=0;s<p.row_NS;s++) {
205     EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=m_t/diagS;
206     }
207     #else /* row-sum lumping */
208     for (s=0;s<p.row_NS;s++) {
209     rtmp=0;
210     for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
211     EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp*D_p[0];
212     }
213     #endif
214     for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
215     Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
216 gross 1204 } /* end color check */
217     } /* end element loop */
218     } /* end color loop */
219     }
220 gross 798 } else {
221 gross 1204 if (expandedD) {
222     #ifndef PASO_MPI
223     for (color=elements->minColor;color<=elements->maxColor;color++) {
224     /* open loop over all elements: */
225     #pragma omp for private(e) schedule(static)
226     for(e=0;e<elements->numElements;e++){
227     if (elements->Color[e]==color) {
228     #else
229     {
230     for(e=0;e<elements->numElements;e++) {
231     {
232     #endif
233     Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
234     memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
235     D_p=getSampleData(D,e);
236 btully 1220 #ifdef NEW_LUMPING /* HRZ lumping */
237     /*
238     * Number of PDEs: Multiple
239     * D_p varies over element: True
240     */
241     for (k=0;k<p.numEqu;k++) {
242     m_t=0; /* mass of the element: m_t */
243     for (q=0;q<p.numQuad;q++) {
244     m_t+=Vol[q]*D_p[INDEX2(k,q,p.numEqu)];
245     }
246     diagS=0; /* diagonal sum: S */
247     for (s=0;s<p.row_NS;s++) {
248     rtmp=0;
249     for (q=0;q<p.numQuad;q++) {
250     rtmp+=Vol[q]*D_p[INDEX2(k,q,p.numEqu)]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
251     }
252     EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
253     diagS+=EM_lumpedMat[INDEX2(k,s,p.numEqu)];
254     }
255     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
256     for (s=0;s<p.row_NS;s++) {
257     EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=m_t/diagS;
258     }
259     }
260     #else /* row-sum lumping */
261 gross 1204 for (s=0;s<p.row_NS;s++) {
262     for (k=0;k<p.numEqu;k++) {
263     rtmp=0.;
264     for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[INDEX2(k,q,p.numEqu)];
265     EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
266     }
267     }
268 btully 1220 #endif
269 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)]];
270     Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
271     } /* end color check */
272     } /* end element loop */
273     } /* end color loop */
274     } else {
275     #ifndef PASO_MPI
276     for (color=elements->minColor;color<=elements->maxColor;color++) {
277     /* open loop over all elements: */
278     #pragma omp for private(e) schedule(static)
279     for(e=0;e<elements->numElements;e++){
280     if (elements->Color[e]==color) {
281     #else
282     {
283     for(e=0;e<elements->numElements;e++) {
284     {
285     #endif
286     Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
287     memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
288     D_p=getSampleData(D,e);
289 btully 1220 #ifdef NEW_LUMPING /* HRZ lumping */
290     /*
291     * Number of PDEs: Multiple
292     * D_p varies over element: False
293     */
294     for (k=0;k<p.numEqu;k++) {
295     m_t=0; /* mass of the element: m_t */
296     for (q=0;q<p.numQuad;q++) {
297     m_t+=Vol[q]*D_p[k];
298     }
299     diagS=0; /* diagonal sum: S */
300     for (s=0;s<p.row_NS;s++) {
301     rtmp=0;
302     for (q=0;q<p.numQuad;q++) {
303     rtmp+=Vol[q]*D_p[k]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
304     }
305     EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
306     diagS+=EM_lumpedMat[INDEX2(k,s,p.numEqu)];
307     }
308     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
309     for (s=0;s<p.row_NS;s++) {
310     EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=m_t/diagS;
311     }
312     }
313     #else /* row-sum lumping */
314 gross 1204 for (s=0;s<p.row_NS;s++) {
315     rtmp=0;
316     for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
317     for (k=0;k<p.numEqu;k++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp*D_p[k];
318     }
319 btully 1220 #endif
320 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)]];
321     Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
322     } /* end color check */
323     } /* end element loop */
324     } /* end color loop */
325     }
326 gross 798 }
327 gross 1204 } /* end of pointer check */
328     THREAD_MEMFREE(EM_lumpedMat);
329     THREAD_MEMFREE(row_index);
330     } /* end parallel region */
331 jgs 82 }
332 gross 1204 #ifdef Finley_TRACE
333     printf("timing: assemblage lumped PDE: %.4e sec\n",Finley_timer()-time0);
334     #endif
335 jgs 82 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26