/[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 1220 - (hide annotations)
Thu Aug 2 04:46:20 2007 UTC (12 years ago) by btully
File MIME type: text/plain
File size: 15898 byte(s)
Included HRZ lumping scheme for the LUMPED solver.

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 btully 1220 #define NEW_LUMPING /* */
37    
38 jgs 82 /**************************************************************/
39    
40 gross 1204 void Finley_Assemble_LumpedSystem(Finley_NodeFile* nodes,Finley_ElementFile* elements, escriptDataC* lumpedMat, escriptDataC* D)
41     {
42 jgs 82
43 gross 1204 bool_t reducedIntegrationOrder=FALSE, expandedD;
44 jgs 150 char error_msg[LenErrorMsg_MAX];
45 gross 798 Assemble_Parameters p;
46 jgs 82 double time0;
47 gross 1204 dim_t dimensions[ESCRIPT_MAX_DATA_RANK], k, e, len_EM_lumpedMat, q, s;
48 gross 1028 type_t funcspace;
49 gross 1204 index_t color,*row_index=NULL;
50     double *S=NULL, *EM_lumpedMat=NULL, *Vol=NULL, *D_p=NULL, *lumpedMat_p=NULL;
51 btully 1220 register double rtmp, m_t, diagS;
52 gross 1204 size_t len_EM_lumpedMat_size;
53 gross 798
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 jgs 82 if (! numSamplesEqual(D,p.numQuad,elements->numElements) ) {
84 gross 1204 sprintf(error_msg,"Assemble_LumpedSystem: sample points of coefficient D don't match (%d,%d)",p.numQuad,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 gross 1204 lumpedMat_p=getSampleData(lumpedMat,0);
109     len_EM_lumpedMat=p.row_NN*p.numEqu;
110     len_EM_lumpedMat_size=len_EM_lumpedMat*sizeof(double);
111     expandedD=isExpanded(D);
112     S=p.row_jac->ReferenceElement->S;
113    
114 btully 1220 #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp, m_t, diagS)
115 gross 1204 {
116     EM_lumpedMat=THREAD_MEMALLOC(len_EM_lumpedMat,double);
117     row_index=THREAD_MEMALLOC(p.row_NN,index_t);
118     if ( !Finley_checkPtr(EM_lumpedMat) && !Finley_checkPtr(row_index) ) {
119     if (p.numEqu == 1) {
120     if (expandedD) {
121     #ifndef PASO_MPI
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     #else
128     {
129     for(e=0;e<elements->numElements;e++) {
130     {
131     #endif
132     Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
133     memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
134     D_p=getSampleData(D,e);
135 btully 1220 #ifdef NEW_LUMPING /* HRZ lumping */
136     /*
137     * Number of PDEs: 1
138     * D_p varies over element: True
139     */
140     m_t=0; /* mass of the element: m_t */
141     for (q=0;q<p.numQuad;q++) {
142     m_t+=Vol[q]*D_p[q];
143     }
144     diagS=0; /* diagonal sum: S */
145 gross 1204 for (s=0;s<p.row_NS;s++) {
146     rtmp=0;
147 btully 1220 for (q=0;q<p.numQuad;q++) {
148     rtmp+=Vol[q]*D_p[q]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
149     }
150     EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
151     diagS+=EM_lumpedMat[INDEX2(0,s,p.numEqu)];
152     }
153     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
154     for (s=0;s<p.row_NS;s++) {
155     EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=m_t/diagS;
156     }
157     #else /* row-sum lumping */
158     for (s=0;s<p.row_NS;s++) {
159     rtmp=0;
160 gross 1204 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[q];
161     EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
162     }
163 btully 1220 #endif
164 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)]];
165     Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
166     } /* end color check */
167     } /* end element loop */
168     } /* end color loop */
169     } else {
170     #ifndef PASO_MPI
171     for (color=elements->minColor;color<=elements->maxColor;color++) {
172     /* open loop over all elements: */
173     #pragma omp for private(e) schedule(static)
174     for(e=0;e<elements->numElements;e++){
175     if (elements->Color[e]==color) {
176     #else
177     {
178     for(e=0;e<elements->numElements;e++) {
179     {
180     #endif
181 btully 1220 Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
182     memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
183     D_p=getSampleData(D,e);
184     #ifdef NEW_LUMPING /* HRZ lumping */
185     /*
186     * Number of PDEs: 1
187     * D_p varies over element: False
188     */
189     m_t=0; /* mass of the element: m_t */
190     for (q=0;q<p.numQuad;q++) {
191     m_t+=Vol[q]*D_p[0];
192     }
193     diagS=0; /* diagonal sum: S */
194     for (s=0;s<p.row_NS;s++) {
195     rtmp=0;
196     for (q=0;q<p.numQuad;q++) {
197     rtmp+=Vol[q]*D_p[0]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
198     }
199     EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
200     diagS+=EM_lumpedMat[INDEX2(0,s,p.numEqu)];
201     }
202     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
203     for (s=0;s<p.row_NS;s++) {
204     EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=m_t/diagS;
205     }
206     #else /* row-sum lumping */
207     for (s=0;s<p.row_NS;s++) {
208     rtmp=0;
209     for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
210     EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp*D_p[0];
211     }
212     #endif
213     for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
214     Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
215 gross 1204 } /* end color check */
216     } /* end element loop */
217     } /* end color loop */
218     }
219 gross 798 } else {
220 gross 1204 if (expandedD) {
221     #ifndef PASO_MPI
222     for (color=elements->minColor;color<=elements->maxColor;color++) {
223     /* open loop over all elements: */
224     #pragma omp for private(e) schedule(static)
225     for(e=0;e<elements->numElements;e++){
226     if (elements->Color[e]==color) {
227     #else
228     {
229     for(e=0;e<elements->numElements;e++) {
230     {
231     #endif
232     Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
233     memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
234     D_p=getSampleData(D,e);
235 btully 1220 #ifdef NEW_LUMPING /* HRZ lumping */
236     /*
237     * Number of PDEs: Multiple
238     * D_p varies over element: True
239     */
240     for (k=0;k<p.numEqu;k++) {
241     m_t=0; /* mass of the element: m_t */
242     for (q=0;q<p.numQuad;q++) {
243     m_t+=Vol[q]*D_p[INDEX2(k,q,p.numEqu)];
244     }
245     diagS=0; /* diagonal sum: S */
246     for (s=0;s<p.row_NS;s++) {
247     rtmp=0;
248     for (q=0;q<p.numQuad;q++) {
249     rtmp+=Vol[q]*D_p[INDEX2(k,q,p.numEqu)]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
250     }
251     EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
252     diagS+=EM_lumpedMat[INDEX2(k,s,p.numEqu)];
253     }
254     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
255     for (s=0;s<p.row_NS;s++) {
256     EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=m_t/diagS;
257     }
258     }
259     #else /* row-sum lumping */
260 gross 1204 for (s=0;s<p.row_NS;s++) {
261     for (k=0;k<p.numEqu;k++) {
262     rtmp=0.;
263     for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[INDEX2(k,q,p.numEqu)];
264     EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
265     }
266     }
267 btully 1220 #endif
268 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)]];
269     Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
270     } /* end color check */
271     } /* end element loop */
272     } /* end color loop */
273     } else {
274     #ifndef PASO_MPI
275     for (color=elements->minColor;color<=elements->maxColor;color++) {
276     /* open loop over all elements: */
277     #pragma omp for private(e) schedule(static)
278     for(e=0;e<elements->numElements;e++){
279     if (elements->Color[e]==color) {
280     #else
281     {
282     for(e=0;e<elements->numElements;e++) {
283     {
284     #endif
285     Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
286     memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
287     D_p=getSampleData(D,e);
288 btully 1220 #ifdef NEW_LUMPING /* HRZ lumping */
289     /*
290     * Number of PDEs: Multiple
291     * D_p varies over element: False
292     */
293     for (k=0;k<p.numEqu;k++) {
294     m_t=0; /* mass of the element: m_t */
295     for (q=0;q<p.numQuad;q++) {
296     m_t+=Vol[q]*D_p[k];
297     }
298     diagS=0; /* diagonal sum: S */
299     for (s=0;s<p.row_NS;s++) {
300     rtmp=0;
301     for (q=0;q<p.numQuad;q++) {
302     rtmp+=Vol[q]*D_p[k]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
303     }
304     EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
305     diagS+=EM_lumpedMat[INDEX2(k,s,p.numEqu)];
306     }
307     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
308     for (s=0;s<p.row_NS;s++) {
309     EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=m_t/diagS;
310     }
311     }
312     #else /* row-sum lumping */
313 gross 1204 for (s=0;s<p.row_NS;s++) {
314     rtmp=0;
315     for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
316     for (k=0;k<p.numEqu;k++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp*D_p[k];
317     }
318 btully 1220 #endif
319 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)]];
320     Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
321     } /* end color check */
322     } /* end element loop */
323     } /* end color loop */
324     }
325 gross 798 }
326 gross 1204 } /* end of pointer check */
327     THREAD_MEMFREE(EM_lumpedMat);
328     THREAD_MEMFREE(row_index);
329     } /* end parallel region */
330 jgs 82 }
331 gross 1204 #ifdef Finley_TRACE
332     printf("timing: assemblage lumped PDE: %.4e sec\n",Finley_timer()-time0);
333     #endif
334 jgs 82 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26