/[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 1650 - (show annotations)
Tue Jul 15 08:58:22 2008 UTC (11 years, 1 month ago) by phornby
File MIME type: text/plain
File size: 14975 byte(s)
more care with m_t and diagS needed in the pragmas.
1
2 /* $Id$ */
3
4 /*******************************************************
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 /**************************************************************/
17
18 /* assembles the mass matrix in lumped form */
19
20 /* The coefficient D has to be defined on the integration points or not present. */
21
22 /* lumpedMat has to be initialized before the routine is called. */
23
24 /**************************************************************/
25
26 #include "Assemble.h"
27 #include "Util.h"
28 #ifdef _OPENMP
29 #include <omp.h>
30 #endif
31
32
33 /* Disabled until the tests pass */
34 /* #define NEW_LUMPING */ /* */
35
36 /**************************************************************/
37
38 void Finley_Assemble_LumpedSystem(Finley_NodeFile* nodes,Finley_ElementFile* elements, escriptDataC* lumpedMat, escriptDataC* D)
39 {
40
41 bool_t reducedIntegrationOrder=FALSE, expandedD;
42 char error_msg[LenErrorMsg_MAX];
43 Assemble_Parameters p;
44 dim_t dimensions[ESCRIPT_MAX_DATA_RANK], k, e, len_EM_lumpedMat, q, s;
45 type_t funcspace;
46 index_t color,*row_index=NULL;
47 double *S=NULL, *EM_lumpedMat=NULL, *Vol=NULL, *D_p=NULL, *lumpedMat_p=NULL;
48 register double rtmp;
49 size_t len_EM_lumpedMat_size;
50
51 #ifdef NEW_LUMPING
52 register double m_t, diagS;
53 #endif
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)
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 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 Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
128 memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
129 D_p=getSampleData(D,e);
130 #ifdef NEW_LUMPING /* HRZ lumping */
131 /*
132 * Number of PDEs: 1
133 * D_p varies over element: True
134 */
135 #pragma omp parallel private(m_t)
136 m_t=0; /* mass of the element: m_t */
137 for (q=0;q<p.numQuad;q++) {
138 m_t+=Vol[q]*D_p[q];
139 }
140 #pragma omp parallel private(diagS)
141 diagS=0; /* diagonal sum: S */
142 for (s=0;s<p.row_NS;s++) {
143 rtmp=0;
144 for (q=0;q<p.numQuad;q++) {
145 rtmp+=Vol[q]*D_p[q]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
146 }
147 EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
148 diagS+=EM_lumpedMat[INDEX2(0,s,p.numEqu)];
149 }
150 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
151 for (s=0;s<p.row_NS;s++) {
152 EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=m_t/diagS;
153 }
154 #else /* row-sum lumping */
155 for (s=0;s<p.row_NS;s++) {
156 rtmp=0;
157 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[q];
158 EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
159 }
160 #endif
161 for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
162 Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
163 } /* end color check */
164 } /* end element loop */
165 } /* end color loop */
166 } else {
167 for (color=elements->minColor;color<=elements->maxColor;color++) {
168 /* open loop over all elements: */
169 #pragma omp for private(e) schedule(static)
170 for(e=0;e<elements->numElements;e++){
171 if (elements->Color[e]==color) {
172 Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
173 memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
174 D_p=getSampleData(D,e);
175 #ifdef NEW_LUMPING /* HRZ lumping */
176 /*
177 * Number of PDEs: 1
178 * D_p varies over element: False
179 */
180 m_t=0; /* mass of the element: m_t */
181 for (q=0;q<p.numQuad;q++) {
182 m_t+=Vol[q]*D_p[0];
183 }
184 diagS=0; /* diagonal sum: S */
185 for (s=0;s<p.row_NS;s++) {
186 rtmp=0;
187 for (q=0;q<p.numQuad;q++) {
188 rtmp+=Vol[q]*D_p[0]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
189 }
190 EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
191 diagS+=EM_lumpedMat[INDEX2(0,s,p.numEqu)];
192 }
193 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
194 for (s=0;s<p.row_NS;s++) {
195 EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=m_t/diagS;
196 }
197 #else /* row-sum lumping */
198 for (s=0;s<p.row_NS;s++) {
199 rtmp=0;
200 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
201 EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp*D_p[0];
202 }
203 #endif
204 for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
205 Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
206 } /* end color check */
207 } /* end element loop */
208 } /* end color loop */
209 }
210 } else {
211 if (expandedD) {
212 for (color=elements->minColor;color<=elements->maxColor;color++) {
213 /* open loop over all elements: */
214 #pragma omp for private(e) schedule(static)
215 for(e=0;e<elements->numElements;e++){
216 if (elements->Color[e]==color) {
217 Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
218 memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
219 D_p=getSampleData(D,e);
220 #ifdef NEW_LUMPING /* HRZ lumping */
221 /*
222 * Number of PDEs: Multiple
223 * D_p varies over element: True
224 */
225 for (k=0;k<p.numEqu;k++) {
226 m_t=0; /* mass of the element: m_t */
227 for (q=0;q<p.numQuad;q++) {
228 m_t+=Vol[q]*D_p[INDEX2(k,q,p.numEqu)];
229 }
230 diagS=0; /* diagonal sum: S */
231 for (s=0;s<p.row_NS;s++) {
232 rtmp=0;
233 for (q=0;q<p.numQuad;q++) {
234 rtmp+=Vol[q]*D_p[INDEX2(k,q,p.numEqu)]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
235 }
236 EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
237 diagS+=EM_lumpedMat[INDEX2(k,s,p.numEqu)];
238 }
239 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
240 for (s=0;s<p.row_NS;s++) {
241 EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=m_t/diagS;
242 }
243 }
244 #else /* row-sum lumping */
245 for (s=0;s<p.row_NS;s++) {
246 for (k=0;k<p.numEqu;k++) {
247 rtmp=0.;
248 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[INDEX2(k,q,p.numEqu)];
249 EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
250 }
251 }
252 #endif
253 for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
254 Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
255 } /* end color check */
256 } /* end element loop */
257 } /* end color loop */
258 } else {
259 /* open loop over all elements: */
260 for (color=elements->minColor;color<=elements->maxColor;color++) {
261 #pragma omp for private(e) schedule(static)
262 for(e=0;e<elements->numElements;e++){
263 if (elements->Color[e]==color) {
264 Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
265 memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
266 D_p=getSampleData(D,e);
267 #ifdef NEW_LUMPING /* HRZ lumping */
268 /*
269 * Number of PDEs: Multiple
270 * D_p varies over element: False
271 */
272 for (k=0;k<p.numEqu;k++) {
273 m_t=0; /* mass of the element: m_t */
274 for (q=0;q<p.numQuad;q++) {
275 m_t+=Vol[q]*D_p[k];
276 }
277 diagS=0; /* diagonal sum: S */
278 for (s=0;s<p.row_NS;s++) {
279 rtmp=0;
280 for (q=0;q<p.numQuad;q++) {
281 rtmp+=Vol[q]*D_p[k]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
282 }
283 EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
284 diagS+=EM_lumpedMat[INDEX2(k,s,p.numEqu)];
285 }
286 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
287 for (s=0;s<p.row_NS;s++) {
288 EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=m_t/diagS;
289 }
290 }
291 #else /* row-sum lumping */
292 for (s=0;s<p.row_NS;s++) {
293 rtmp=0;
294 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
295 for (k=0;k<p.numEqu;k++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp*D_p[k];
296 }
297 #endif
298 for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
299 Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
300 } /* end color check */
301 } /* end element loop */
302 } /* end color loop */
303 }
304 }
305 } /* end of pointer check */
306 THREAD_MEMFREE(EM_lumpedMat);
307 THREAD_MEMFREE(row_index);
308 } /* end parallel region */
309 }
310 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26