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

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 #define NEW_LUMPING /* */
37
38 /**************************************************************/
39
40 void Finley_Assemble_LumpedSystem(Finley_NodeFile* nodes,Finley_ElementFile* elements, escriptDataC* lumpedMat, escriptDataC* D)
41 {
42
43 bool_t reducedIntegrationOrder=FALSE, expandedD;
44 char error_msg[LenErrorMsg_MAX];
45 Assemble_Parameters p;
46 double time0;
47 dim_t dimensions[ESCRIPT_MAX_DATA_RANK], k, e, len_EM_lumpedMat, q, s;
48 type_t funcspace;
49 index_t color,*row_index=NULL;
50 double *S=NULL, *EM_lumpedMat=NULL, *Vol=NULL, *D_p=NULL, *lumpedMat_p=NULL;
51 register double rtmp, m_t, diagS;
52 size_t len_EM_lumpedMat_size;
53
54 Finley_resetError();
55
56 if (nodes==NULL || elements==NULL) return;
57 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 }
62 funcspace=getFunctionSpaceType(D);
63 /* check if all function spaces are the same */
64 if (funcspace==FINLEY_ELEMENTS) {
65 reducedIntegrationOrder=FALSE;
66 } else if (funcspace==FINLEY_FACE_ELEMENTS) {
67 reducedIntegrationOrder=FALSE;
68 } else if (funcspace==FINLEY_REDUCED_ELEMENTS) {
69 reducedIntegrationOrder=TRUE;
70 } else if (funcspace==FINLEY_REDUCED_FACE_ELEMENTS) {
71 reducedIntegrationOrder=TRUE;
72 } else {
73 Finley_setError(TYPE_ERROR,"Assemble_LumpedSystem: assemblage failed because of illegal function space.");
74 }
75 if (! Finley_noError()) return;
76
77 /* set all parameters in p*/
78 Assemble_getAssembleParameters(nodes,elements,NULL,lumpedMat, reducedIntegrationOrder, &p);
79 if (! Finley_noError()) return;
80
81 /* check if all function spaces are the same */
82
83 if (! numSamplesEqual(D,p.numQuad,elements->numElements) ) {
84 sprintf(error_msg,"Assemble_LumpedSystem: sample points of coefficient D don't match (%d,%d)",p.numQuad,elements->numElements);
85 Finley_setError(TYPE_ERROR,error_msg);
86 }
87
88 /* check the dimensions: */
89
90 if (p.numEqu==1) {
91 if (!isEmpty(D)) {
92 if (!isDataPointShapeEqual(D,0,dimensions)) {
93 Finley_setError(TYPE_ERROR,"Assemble_LumpedSystem: coefficient D, rank 0 expected.");
94 }
95
96 }
97 } else {
98 if (!isEmpty(D)) {
99 dimensions[0]=p.numEqu;
100 if (!isDataPointShapeEqual(D,1,dimensions)) {
101 sprintf(error_msg,"Assemble_LumpedSystem: coefficient D, expected shape (%d,)",dimensions[0]);
102 Finley_setError(TYPE_ERROR,error_msg);
103 }
104 }
105 }
106
107 if (Finley_noError()) {
108 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 #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp, m_t, diagS)
115 {
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 #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 for (s=0;s<p.row_NS;s++) {
146 rtmp=0;
147 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 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 #endif
164 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 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 } /* end color check */
216 } /* end element loop */
217 } /* end color loop */
218 }
219 } else {
220 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 #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 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 #endif
268 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 #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 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 #endif
319 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 }
326 } /* end of pointer check */
327 THREAD_MEMFREE(EM_lumpedMat);
328 THREAD_MEMFREE(row_index);
329 } /* end parallel region */
330 }
331 #ifdef Finley_TRACE
332 printf("timing: assemblage lumped PDE: %.4e sec\n",Finley_timer()-time0);
333 #endif
334 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26