/[escript]/branches/arrexp_2137_win_merge/finley/src/Assemble_LumpedSystem.c
ViewVC logotype

Contents of /branches/arrexp_2137_win_merge/finley/src/Assemble_LumpedSystem.c

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26