/[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 1312 - (hide annotations)
Mon Sep 24 06:18:44 2007 UTC (11 years, 3 months ago) by ksteube
File MIME type: text/plain
File size: 14943 byte(s)
The MPI branch is hereby closed. All future work should be in trunk.

Previously in revision 1295 I merged the latest changes to trunk into trunk-mpi-branch.
In this revision I copied all files from trunk-mpi-branch over the corresponding
trunk files. I did not use 'svn merge', it was a copy.

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 jgs 82 double time0;
45 gross 1204 dim_t dimensions[ESCRIPT_MAX_DATA_RANK], k, e, len_EM_lumpedMat, q, s;
46 gross 1028 type_t funcspace;
47 gross 1204 index_t color,*row_index=NULL;
48     double *S=NULL, *EM_lumpedMat=NULL, *Vol=NULL, *D_p=NULL, *lumpedMat_p=NULL;
49 btully 1220 register double rtmp, m_t, diagS;
50 gross 1204 size_t len_EM_lumpedMat_size;
51 gross 798
52 jgs 150 Finley_resetError();
53 jgs 82
54     if (nodes==NULL || elements==NULL) return;
55 gross 1204 if (isEmpty(lumpedMat) || isEmpty(D)) return;
56     if (isEmpty(lumpedMat) && !isEmpty(D)) {
57     Finley_setError(TYPE_ERROR,"Assemble_LumpedSystem: coefficients are non-zero but no lumped matrix is given.");
58     return;
59 gross 798 }
60 gross 1204 funcspace=getFunctionSpaceType(D);
61 jgs 82 /* check if all function spaces are the same */
62 gross 798 if (funcspace==FINLEY_ELEMENTS) {
63     reducedIntegrationOrder=FALSE;
64     } else if (funcspace==FINLEY_FACE_ELEMENTS) {
65     reducedIntegrationOrder=FALSE;
66 gross 1072 } else if (funcspace==FINLEY_REDUCED_ELEMENTS) {
67     reducedIntegrationOrder=TRUE;
68     } else if (funcspace==FINLEY_REDUCED_FACE_ELEMENTS) {
69     reducedIntegrationOrder=TRUE;
70 gross 798 } else {
71 gross 1204 Finley_setError(TYPE_ERROR,"Assemble_LumpedSystem: assemblage failed because of illegal function space.");
72 gross 798 }
73     if (! Finley_noError()) return;
74 jgs 82
75 gross 798 /* set all parameters in p*/
76 gross 1204 Assemble_getAssembleParameters(nodes,elements,NULL,lumpedMat, reducedIntegrationOrder, &p);
77 gross 798 if (! Finley_noError()) return;
78    
79     /* check if all function spaces are the same */
80    
81 jgs 82 if (! numSamplesEqual(D,p.numQuad,elements->numElements) ) {
82 gross 1204 sprintf(error_msg,"Assemble_LumpedSystem: sample points of coefficient D don't match (%d,%d)",p.numQuad,elements->numElements);
83 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
84 jgs 82 }
85    
86     /* check the dimensions: */
87    
88 gross 1204 if (p.numEqu==1) {
89 jgs 82 if (!isEmpty(D)) {
90     if (!isDataPointShapeEqual(D,0,dimensions)) {
91 gross 1204 Finley_setError(TYPE_ERROR,"Assemble_LumpedSystem: coefficient D, rank 0 expected.");
92 jgs 82 }
93 gross 1204
94 jgs 82 }
95     } else {
96     if (!isEmpty(D)) {
97     dimensions[0]=p.numEqu;
98 gross 1204 if (!isDataPointShapeEqual(D,1,dimensions)) {
99     sprintf(error_msg,"Assemble_LumpedSystem: coefficient D, expected shape (%d,)",dimensions[0]);
100 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
101 jgs 82 }
102     }
103     }
104 gross 1204
105 jgs 150 if (Finley_noError()) {
106 gross 1204 lumpedMat_p=getSampleData(lumpedMat,0);
107     len_EM_lumpedMat=p.row_NN*p.numEqu;
108     len_EM_lumpedMat_size=len_EM_lumpedMat*sizeof(double);
109     expandedD=isExpanded(D);
110     S=p.row_jac->ReferenceElement->S;
111    
112 btully 1220 #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp, m_t, diagS)
113 gross 1204 {
114     EM_lumpedMat=THREAD_MEMALLOC(len_EM_lumpedMat,double);
115     row_index=THREAD_MEMALLOC(p.row_NN,index_t);
116     if ( !Finley_checkPtr(EM_lumpedMat) && !Finley_checkPtr(row_index) ) {
117     if (p.numEqu == 1) {
118     if (expandedD) {
119     for (color=elements->minColor;color<=elements->maxColor;color++) {
120     /* open loop over all elements: */
121     #pragma omp for private(e) schedule(static)
122     for(e=0;e<elements->numElements;e++){
123     if (elements->Color[e]==color) {
124     Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
125     memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
126     D_p=getSampleData(D,e);
127 btully 1220 #ifdef NEW_LUMPING /* HRZ lumping */
128     /*
129     * Number of PDEs: 1
130     * D_p varies over element: True
131     */
132     m_t=0; /* mass of the element: m_t */
133     for (q=0;q<p.numQuad;q++) {
134     m_t+=Vol[q]*D_p[q];
135     }
136     diagS=0; /* diagonal sum: S */
137 gross 1204 for (s=0;s<p.row_NS;s++) {
138     rtmp=0;
139 btully 1220 for (q=0;q<p.numQuad;q++) {
140     rtmp+=Vol[q]*D_p[q]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
141     }
142     EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
143     diagS+=EM_lumpedMat[INDEX2(0,s,p.numEqu)];
144     }
145     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
146     for (s=0;s<p.row_NS;s++) {
147     EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=m_t/diagS;
148     }
149     #else /* row-sum lumping */
150     for (s=0;s<p.row_NS;s++) {
151     rtmp=0;
152 gross 1204 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[q];
153     EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
154     }
155 btully 1220 #endif
156 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)]];
157     Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
158     } /* end color check */
159     } /* end element loop */
160     } /* end color loop */
161     } else {
162     for (color=elements->minColor;color<=elements->maxColor;color++) {
163     /* open loop over all elements: */
164     #pragma omp for private(e) schedule(static)
165     for(e=0;e<elements->numElements;e++){
166     if (elements->Color[e]==color) {
167 btully 1220 Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
168     memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
169     D_p=getSampleData(D,e);
170     #ifdef NEW_LUMPING /* HRZ lumping */
171     /*
172     * Number of PDEs: 1
173     * D_p varies over element: False
174     */
175     m_t=0; /* mass of the element: m_t */
176     for (q=0;q<p.numQuad;q++) {
177     m_t+=Vol[q]*D_p[0];
178     }
179     diagS=0; /* diagonal sum: S */
180     for (s=0;s<p.row_NS;s++) {
181     rtmp=0;
182     for (q=0;q<p.numQuad;q++) {
183     rtmp+=Vol[q]*D_p[0]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
184     }
185     EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
186     diagS+=EM_lumpedMat[INDEX2(0,s,p.numEqu)];
187     }
188     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
189     for (s=0;s<p.row_NS;s++) {
190     EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=m_t/diagS;
191     }
192     #else /* row-sum lumping */
193     for (s=0;s<p.row_NS;s++) {
194     rtmp=0;
195     for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
196     EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp*D_p[0];
197     }
198     #endif
199     for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
200     Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
201 gross 1204 } /* end color check */
202     } /* end element loop */
203     } /* end color loop */
204     }
205 gross 798 } else {
206 gross 1204 if (expandedD) {
207     for (color=elements->minColor;color<=elements->maxColor;color++) {
208     /* open loop over all elements: */
209     #pragma omp for private(e) schedule(static)
210     for(e=0;e<elements->numElements;e++){
211     if (elements->Color[e]==color) {
212     Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
213     memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
214     D_p=getSampleData(D,e);
215 btully 1220 #ifdef NEW_LUMPING /* HRZ lumping */
216     /*
217     * Number of PDEs: Multiple
218     * D_p varies over element: True
219     */
220     for (k=0;k<p.numEqu;k++) {
221     m_t=0; /* mass of the element: m_t */
222     for (q=0;q<p.numQuad;q++) {
223     m_t+=Vol[q]*D_p[INDEX2(k,q,p.numEqu)];
224     }
225     diagS=0; /* diagonal sum: S */
226     for (s=0;s<p.row_NS;s++) {
227     rtmp=0;
228     for (q=0;q<p.numQuad;q++) {
229     rtmp+=Vol[q]*D_p[INDEX2(k,q,p.numEqu)]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
230     }
231     EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
232     diagS+=EM_lumpedMat[INDEX2(k,s,p.numEqu)];
233     }
234     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
235     for (s=0;s<p.row_NS;s++) {
236     EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=m_t/diagS;
237     }
238     }
239     #else /* row-sum lumping */
240 gross 1204 for (s=0;s<p.row_NS;s++) {
241     for (k=0;k<p.numEqu;k++) {
242     rtmp=0.;
243     for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[INDEX2(k,q,p.numEqu)];
244     EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
245     }
246     }
247 btully 1220 #endif
248 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)]];
249     Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
250     } /* end color check */
251     } /* end element loop */
252     } /* end color loop */
253     } else {
254 ksteube 1312 /* open loop over all elements: */
255 gross 1204 for (color=elements->minColor;color<=elements->maxColor;color++) {
256     #pragma omp for private(e) schedule(static)
257     for(e=0;e<elements->numElements;e++){
258     if (elements->Color[e]==color) {
259     Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
260     memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
261     D_p=getSampleData(D,e);
262 btully 1220 #ifdef NEW_LUMPING /* HRZ lumping */
263     /*
264     * Number of PDEs: Multiple
265     * D_p varies over element: False
266     */
267     for (k=0;k<p.numEqu;k++) {
268     m_t=0; /* mass of the element: m_t */
269     for (q=0;q<p.numQuad;q++) {
270     m_t+=Vol[q]*D_p[k];
271     }
272     diagS=0; /* diagonal sum: S */
273     for (s=0;s<p.row_NS;s++) {
274     rtmp=0;
275     for (q=0;q<p.numQuad;q++) {
276     rtmp+=Vol[q]*D_p[k]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
277     }
278     EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
279     diagS+=EM_lumpedMat[INDEX2(k,s,p.numEqu)];
280     }
281     /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
282     for (s=0;s<p.row_NS;s++) {
283     EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=m_t/diagS;
284     }
285     }
286     #else /* row-sum lumping */
287 gross 1204 for (s=0;s<p.row_NS;s++) {
288     rtmp=0;
289     for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
290     for (k=0;k<p.numEqu;k++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp*D_p[k];
291     }
292 btully 1220 #endif
293 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)]];
294     Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
295     } /* end color check */
296     } /* end element loop */
297     } /* end color loop */
298     }
299 gross 798 }
300 gross 1204 } /* end of pointer check */
301     THREAD_MEMFREE(EM_lumpedMat);
302     THREAD_MEMFREE(row_index);
303     } /* end parallel region */
304 jgs 82 }
305 gross 1204 #ifdef Finley_TRACE
306     printf("timing: assemblage lumped PDE: %.4e sec\n",Finley_timer()-time0);
307     #endif
308 jgs 82 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26