/[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 1811 - (hide annotations)
Thu Sep 25 23:11:13 2008 UTC (10 years, 5 months ago) by ksteube
File MIME type: text/plain
File size: 14940 byte(s)
Copyright updated in all files

1 jgs 82
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4     * Copyright (c) 2003-2008 by University of Queensland
5     * 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     double *S=NULL, *EM_lumpedMat=NULL, *Vol=NULL, *D_p=NULL, *lumpedMat_p=NULL;
47 phornby 1628 register double rtmp;
48 gross 1204 size_t len_EM_lumpedMat_size;
49 phornby 1645
50     #ifdef NEW_LUMPING
51 phornby 1628 register double m_t, diagS;
52 phornby 1645 #endif
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 phornby 1650 #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp)
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     for (color=elements->minColor;color<=elements->maxColor;color++) {
122     /* open loop over all elements: */
123     #pragma omp for private(e) schedule(static)
124     for(e=0;e<elements->numElements;e++){
125     if (elements->Color[e]==color) {
126     Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
127     memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
128     D_p=getSampleData(D,e);
129 btully 1220 #ifdef NEW_LUMPING /* HRZ lumping */
130     /*
131     * Number of PDEs: 1
132     * D_p varies over element: True
133     */
134 phornby 1650 #pragma omp parallel private(m_t)
135 btully 1220 m_t=0; /* mass of the element: m_t */
136     for (q=0;q<p.numQuad;q++) {
137     m_t+=Vol[q]*D_p[q];
138     }
139 phornby 1650 #pragma omp parallel private(diagS)
140 btully 1220 diagS=0; /* diagonal sum: S */
141 gross 1204 for (s=0;s<p.row_NS;s++) {
142     rtmp=0;
143 btully 1220 for (q=0;q<p.numQuad;q++) {
144     rtmp+=Vol[q]*D_p[q]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
145     }
146     EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
147     diagS+=EM_lumpedMat[INDEX2(0,s,p.numEqu)];
148     }
149     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
150     for (s=0;s<p.row_NS;s++) {
151     EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=m_t/diagS;
152     }
153     #else /* row-sum lumping */
154     for (s=0;s<p.row_NS;s++) {
155     rtmp=0;
156 gross 1204 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[q];
157     EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
158     }
159 btully 1220 #endif
160 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)]];
161     Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
162     } /* end color check */
163     } /* end element loop */
164     } /* end color loop */
165     } else {
166     for (color=elements->minColor;color<=elements->maxColor;color++) {
167     /* open loop over all elements: */
168     #pragma omp for private(e) schedule(static)
169     for(e=0;e<elements->numElements;e++){
170     if (elements->Color[e]==color) {
171 btully 1220 Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
172     memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
173     D_p=getSampleData(D,e);
174     #ifdef NEW_LUMPING /* HRZ lumping */
175     /*
176     * Number of PDEs: 1
177     * D_p varies over element: False
178     */
179     m_t=0; /* mass of the element: m_t */
180     for (q=0;q<p.numQuad;q++) {
181     m_t+=Vol[q]*D_p[0];
182     }
183     diagS=0; /* diagonal sum: S */
184     for (s=0;s<p.row_NS;s++) {
185     rtmp=0;
186     for (q=0;q<p.numQuad;q++) {
187     rtmp+=Vol[q]*D_p[0]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
188     }
189     EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
190     diagS+=EM_lumpedMat[INDEX2(0,s,p.numEqu)];
191     }
192     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
193     for (s=0;s<p.row_NS;s++) {
194     EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=m_t/diagS;
195     }
196     #else /* row-sum lumping */
197     for (s=0;s<p.row_NS;s++) {
198     rtmp=0;
199     for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
200     EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp*D_p[0];
201     }
202     #endif
203     for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
204     Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
205 gross 1204 } /* end color check */
206     } /* end element loop */
207     } /* end color loop */
208     }
209 gross 798 } else {
210 gross 1204 if (expandedD) {
211     for (color=elements->minColor;color<=elements->maxColor;color++) {
212     /* open loop over all elements: */
213     #pragma omp for private(e) schedule(static)
214     for(e=0;e<elements->numElements;e++){
215     if (elements->Color[e]==color) {
216     Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
217     memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
218     D_p=getSampleData(D,e);
219 btully 1220 #ifdef NEW_LUMPING /* HRZ lumping */
220     /*
221     * Number of PDEs: Multiple
222     * D_p varies over element: True
223     */
224     for (k=0;k<p.numEqu;k++) {
225     m_t=0; /* mass of the element: m_t */
226     for (q=0;q<p.numQuad;q++) {
227     m_t+=Vol[q]*D_p[INDEX2(k,q,p.numEqu)];
228     }
229     diagS=0; /* diagonal sum: S */
230     for (s=0;s<p.row_NS;s++) {
231     rtmp=0;
232     for (q=0;q<p.numQuad;q++) {
233     rtmp+=Vol[q]*D_p[INDEX2(k,q,p.numEqu)]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
234     }
235     EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
236     diagS+=EM_lumpedMat[INDEX2(k,s,p.numEqu)];
237     }
238     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
239     for (s=0;s<p.row_NS;s++) {
240     EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=m_t/diagS;
241     }
242     }
243     #else /* row-sum lumping */
244 gross 1204 for (s=0;s<p.row_NS;s++) {
245     for (k=0;k<p.numEqu;k++) {
246     rtmp=0.;
247     for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[INDEX2(k,q,p.numEqu)];
248     EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
249     }
250     }
251 btully 1220 #endif
252 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)]];
253     Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
254     } /* end color check */
255     } /* end element loop */
256     } /* end color loop */
257     } else {
258 ksteube 1312 /* open loop over all elements: */
259 gross 1204 for (color=elements->minColor;color<=elements->maxColor;color++) {
260     #pragma omp for private(e) schedule(static)
261     for(e=0;e<elements->numElements;e++){
262     if (elements->Color[e]==color) {
263     Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
264     memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
265     D_p=getSampleData(D,e);
266 btully 1220 #ifdef NEW_LUMPING /* HRZ lumping */
267     /*
268     * Number of PDEs: Multiple
269     * D_p varies over element: False
270     */
271     for (k=0;k<p.numEqu;k++) {
272     m_t=0; /* mass of the element: m_t */
273     for (q=0;q<p.numQuad;q++) {
274     m_t+=Vol[q]*D_p[k];
275     }
276     diagS=0; /* diagonal sum: S */
277     for (s=0;s<p.row_NS;s++) {
278     rtmp=0;
279     for (q=0;q<p.numQuad;q++) {
280     rtmp+=Vol[q]*D_p[k]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
281     }
282     EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
283     diagS+=EM_lumpedMat[INDEX2(k,s,p.numEqu)];
284     }
285     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
286     for (s=0;s<p.row_NS;s++) {
287     EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=m_t/diagS;
288     }
289     }
290     #else /* row-sum lumping */
291 gross 1204 for (s=0;s<p.row_NS;s++) {
292     rtmp=0;
293     for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
294     for (k=0;k<p.numEqu;k++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp*D_p[k];
295     }
296 btully 1220 #endif
297 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)]];
298     Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
299     } /* end color check */
300     } /* end element loop */
301     } /* end color loop */
302     }
303 gross 798 }
304 gross 1204 } /* end of pointer check */
305     THREAD_MEMFREE(EM_lumpedMat);
306     THREAD_MEMFREE(row_index);
307     } /* end parallel region */
308 jgs 82 }
309     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26