/[escript]/trunk/finley/src/Assemble_LumpedSystem.c
ViewVC logotype

Contents of /trunk/finley/src/Assemble_LumpedSystem.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1221 - (show 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 /*
2 ************************************************************
3 * Copyright 2006,2007 by ACcENULLNULL MNRF *
4 * *
5 * http://www.access.edu.au *
6 * Primary Business: Queensland, Australia *
7 * Licensed under the Open NULLoftware License version 3.0 *
8 * http://www.opensource.org/licenses/osl-3.0.php *
9 * *
10 ************************************************************
11 */
12
13
14 /**************************************************************/
15
16 /* assembles the mass matrix in lumped form */
17
18 /* The coefficient D has to be defined on the integration points or not present. */
19
20 /* lumpedMat has to be initialized before the routine is called. */
21
22 /**************************************************************/
23
24 /* Author: gross@access.edu.au */
25 /* Version: $Id$ */
26
27 /**************************************************************/
28
29 #include "Assemble.h"
30 #include "Util.h"
31 #ifdef _OPENMP
32 #include <omp.h>
33 #endif
34
35
36 /* Disabled until the tests pass */
37 /* #define NEW_LUMPING */ /* */
38
39 /**************************************************************/
40
41 void Finley_Assemble_LumpedSystem(Finley_NodeFile* nodes,Finley_ElementFile* elements, escriptDataC* lumpedMat, escriptDataC* D)
42 {
43
44 bool_t reducedIntegrationOrder=FALSE, expandedD;
45 char error_msg[LenErrorMsg_MAX];
46 Assemble_Parameters p;
47 double time0;
48 dim_t dimensions[ESCRIPT_MAX_DATA_RANK], k, e, len_EM_lumpedMat, q, s;
49 type_t funcspace;
50 index_t color,*row_index=NULL;
51 double *S=NULL, *EM_lumpedMat=NULL, *Vol=NULL, *D_p=NULL, *lumpedMat_p=NULL;
52 register double rtmp, m_t, diagS;
53 size_t len_EM_lumpedMat_size;
54
55 Finley_resetError();
56
57 if (nodes==NULL || elements==NULL) return;
58 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 }
63 funcspace=getFunctionSpaceType(D);
64 /* check if all function spaces are the same */
65 if (funcspace==FINLEY_ELEMENTS) {
66 reducedIntegrationOrder=FALSE;
67 } else if (funcspace==FINLEY_FACE_ELEMENTS) {
68 reducedIntegrationOrder=FALSE;
69 } else if (funcspace==FINLEY_REDUCED_ELEMENTS) {
70 reducedIntegrationOrder=TRUE;
71 } else if (funcspace==FINLEY_REDUCED_FACE_ELEMENTS) {
72 reducedIntegrationOrder=TRUE;
73 } else {
74 Finley_setError(TYPE_ERROR,"Assemble_LumpedSystem: assemblage failed because of illegal function space.");
75 }
76 if (! Finley_noError()) return;
77
78 /* set all parameters in p*/
79 Assemble_getAssembleParameters(nodes,elements,NULL,lumpedMat, reducedIntegrationOrder, &p);
80 if (! Finley_noError()) return;
81
82 /* check if all function spaces are the same */
83
84 if (! numSamplesEqual(D,p.numQuad,elements->numElements) ) {
85 sprintf(error_msg,"Assemble_LumpedSystem: sample points of coefficient D don't match (%d,%d)",p.numQuad,elements->numElements);
86 Finley_setError(TYPE_ERROR,error_msg);
87 }
88
89 /* check the dimensions: */
90
91 if (p.numEqu==1) {
92 if (!isEmpty(D)) {
93 if (!isDataPointShapeEqual(D,0,dimensions)) {
94 Finley_setError(TYPE_ERROR,"Assemble_LumpedSystem: coefficient D, rank 0 expected.");
95 }
96
97 }
98 } else {
99 if (!isEmpty(D)) {
100 dimensions[0]=p.numEqu;
101 if (!isDataPointShapeEqual(D,1,dimensions)) {
102 sprintf(error_msg,"Assemble_LumpedSystem: coefficient D, expected shape (%d,)",dimensions[0]);
103 Finley_setError(TYPE_ERROR,error_msg);
104 }
105 }
106 }
107
108 if (Finley_noError()) {
109 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 #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp, m_t, diagS)
116 {
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 #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 for (s=0;s<p.row_NS;s++) {
147 rtmp=0;
148 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 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 #endif
165 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 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 } /* end color check */
217 } /* end element loop */
218 } /* end color loop */
219 }
220 } else {
221 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 #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 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 #endif
269 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 #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 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 #endif
320 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 }
327 } /* end of pointer check */
328 THREAD_MEMFREE(EM_lumpedMat);
329 THREAD_MEMFREE(row_index);
330 } /* end parallel region */
331 }
332 #ifdef Finley_TRACE
333 printf("timing: assemblage lumped PDE: %.4e sec\n",Finley_timer()-time0);
334 #endif
335 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26