/[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 2548 - (hide annotations)
Mon Jul 20 06:20:06 2009 UTC (10 years, 1 month ago) by jfenwick
File MIME type: text/plain
File size: 15092 byte(s)
Updating copyright notices
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 1204 dim_t dimensions[ESCRIPT_MAX_DATA_RANK], k, e, len_EM_lumpedMat, q, s;
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 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 jfenwick 2271 void* buffer=allocSampleBuffer(D);
110     requireWrite(lumpedMat);
111     lumpedMat_p=getSampleDataRW(lumpedMat,0);
112 gross 1204 len_EM_lumpedMat=p.row_NN*p.numEqu;
113     len_EM_lumpedMat_size=len_EM_lumpedMat*sizeof(double);
114     expandedD=isExpanded(D);
115     S=p.row_jac->ReferenceElement->S;
116 phornby 1650 #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp)
117 gross 1204 {
118     EM_lumpedMat=THREAD_MEMALLOC(len_EM_lumpedMat,double);
119     row_index=THREAD_MEMALLOC(p.row_NN,index_t);
120     if ( !Finley_checkPtr(EM_lumpedMat) && !Finley_checkPtr(row_index) ) {
121     if (p.numEqu == 1) {
122     if (expandedD) {
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     Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
129     memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
130 jfenwick 2271 D_p=getSampleDataRO(D,e,buffer);
131 btully 1220 #ifdef NEW_LUMPING /* HRZ lumping */
132     /*
133     * Number of PDEs: 1
134     * D_p varies over element: True
135     */
136 phornby 1650 #pragma omp parallel private(m_t)
137 btully 1220 m_t=0; /* mass of the element: m_t */
138     for (q=0;q<p.numQuad;q++) {
139     m_t+=Vol[q]*D_p[q];
140     }
141 phornby 1650 #pragma omp parallel private(diagS)
142 btully 1220 diagS=0; /* diagonal sum: S */
143 gross 1204 for (s=0;s<p.row_NS;s++) {
144     rtmp=0;
145 btully 1220 for (q=0;q<p.numQuad;q++) {
146     rtmp+=Vol[q]*D_p[q]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
147     }
148     EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
149     diagS+=EM_lumpedMat[INDEX2(0,s,p.numEqu)];
150     }
151     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
152     for (s=0;s<p.row_NS;s++) {
153     EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=m_t/diagS;
154     }
155     #else /* row-sum lumping */
156     for (s=0;s<p.row_NS;s++) {
157     rtmp=0;
158 gross 1204 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[q];
159     EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
160     }
161 btully 1220 #endif
162 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)]];
163     Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
164     } /* end color check */
165     } /* end element loop */
166     } /* end color loop */
167     } else {
168     for (color=elements->minColor;color<=elements->maxColor;color++) {
169     /* open loop over all elements: */
170     #pragma omp for private(e) schedule(static)
171     for(e=0;e<elements->numElements;e++){
172     if (elements->Color[e]==color) {
173 btully 1220 Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
174     memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
175 jfenwick 2271 D_p=getSampleDataRO(D,e,buffer);
176 btully 1220 #ifdef NEW_LUMPING /* HRZ lumping */
177     /*
178     * Number of PDEs: 1
179     * D_p varies over element: False
180     */
181     m_t=0; /* mass of the element: m_t */
182     for (q=0;q<p.numQuad;q++) {
183     m_t+=Vol[q]*D_p[0];
184     }
185     diagS=0; /* diagonal sum: S */
186     for (s=0;s<p.row_NS;s++) {
187     rtmp=0;
188     for (q=0;q<p.numQuad;q++) {
189     rtmp+=Vol[q]*D_p[0]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
190     }
191     EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
192     diagS+=EM_lumpedMat[INDEX2(0,s,p.numEqu)];
193     }
194     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
195     for (s=0;s<p.row_NS;s++) {
196     EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=m_t/diagS;
197     }
198     #else /* row-sum lumping */
199     for (s=0;s<p.row_NS;s++) {
200     rtmp=0;
201     for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
202     EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp*D_p[0];
203     }
204     #endif
205     for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
206     Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
207 gross 1204 } /* end color check */
208     } /* end element loop */
209     } /* end color loop */
210     }
211 gross 798 } else {
212 gross 1204 if (expandedD) {
213     for (color=elements->minColor;color<=elements->maxColor;color++) {
214     /* open loop over all elements: */
215     #pragma omp for private(e) schedule(static)
216     for(e=0;e<elements->numElements;e++){
217     if (elements->Color[e]==color) {
218     Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
219     memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
220 jfenwick 2271 D_p=getSampleDataRO(D,e,buffer);
221 btully 1220 #ifdef NEW_LUMPING /* HRZ lumping */
222     /*
223     * Number of PDEs: Multiple
224     * D_p varies over element: True
225     */
226     for (k=0;k<p.numEqu;k++) {
227     m_t=0; /* mass of the element: m_t */
228     for (q=0;q<p.numQuad;q++) {
229     m_t+=Vol[q]*D_p[INDEX2(k,q,p.numEqu)];
230     }
231     diagS=0; /* diagonal sum: S */
232     for (s=0;s<p.row_NS;s++) {
233     rtmp=0;
234     for (q=0;q<p.numQuad;q++) {
235     rtmp+=Vol[q]*D_p[INDEX2(k,q,p.numEqu)]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
236     }
237     EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
238     diagS+=EM_lumpedMat[INDEX2(k,s,p.numEqu)];
239     }
240     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
241     for (s=0;s<p.row_NS;s++) {
242     EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=m_t/diagS;
243     }
244     }
245     #else /* row-sum lumping */
246 gross 1204 for (s=0;s<p.row_NS;s++) {
247     for (k=0;k<p.numEqu;k++) {
248     rtmp=0.;
249     for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[INDEX2(k,q,p.numEqu)];
250     EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
251     }
252     }
253 btully 1220 #endif
254 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)]];
255     Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
256     } /* end color check */
257     } /* end element loop */
258     } /* end color loop */
259     } else {
260 ksteube 1312 /* open loop over all elements: */
261 gross 1204 for (color=elements->minColor;color<=elements->maxColor;color++) {
262     #pragma omp for private(e) schedule(static)
263     for(e=0;e<elements->numElements;e++){
264     if (elements->Color[e]==color) {
265     Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
266     memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
267 jfenwick 2271 D_p=getSampleDataRO(D,e,buffer);
268 btully 1220 #ifdef NEW_LUMPING /* HRZ lumping */
269     /*
270     * Number of PDEs: Multiple
271     * D_p varies over element: False
272     */
273     for (k=0;k<p.numEqu;k++) {
274     m_t=0; /* mass of the element: m_t */
275     for (q=0;q<p.numQuad;q++) {
276     m_t+=Vol[q]*D_p[k];
277     }
278     diagS=0; /* diagonal sum: S */
279     for (s=0;s<p.row_NS;s++) {
280     rtmp=0;
281     for (q=0;q<p.numQuad;q++) {
282     rtmp+=Vol[q]*D_p[k]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
283     }
284     EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
285     diagS+=EM_lumpedMat[INDEX2(k,s,p.numEqu)];
286     }
287     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
288     for (s=0;s<p.row_NS;s++) {
289     EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=m_t/diagS;
290     }
291     }
292     #else /* row-sum lumping */
293 gross 1204 for (s=0;s<p.row_NS;s++) {
294     rtmp=0;
295     for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
296     for (k=0;k<p.numEqu;k++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp*D_p[k];
297     }
298 btully 1220 #endif
299 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)]];
300     Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
301     } /* end color check */
302     } /* end element loop */
303     } /* end color loop */
304     }
305 gross 798 }
306 gross 1204 } /* end of pointer check */
307     THREAD_MEMFREE(EM_lumpedMat);
308     THREAD_MEMFREE(row_index);
309     } /* end parallel region */
310 jfenwick 2271 freeSampleBuffer(buffer);
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