/[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 1384 - (show annotations)
Fri Jan 11 02:29:38 2008 UTC (11 years, 2 months ago) by phornby
Original Path: temp_trunk_copy/finley/src/Assemble_LumpedSystem.c
File MIME type: text/plain
File size: 14943 byte(s)
Make a temp copy of the trunk before checking in the windows changes


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 double time0;
45 dim_t dimensions[ESCRIPT_MAX_DATA_RANK], k, e, len_EM_lumpedMat, q, s;
46 type_t funcspace;
47 index_t color,*row_index=NULL;
48 double *S=NULL, *EM_lumpedMat=NULL, *Vol=NULL, *D_p=NULL, *lumpedMat_p=NULL;
49 register double rtmp, m_t, diagS;
50 size_t len_EM_lumpedMat_size;
51
52 Finley_resetError();
53
54 if (nodes==NULL || elements==NULL) return;
55 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 }
60 funcspace=getFunctionSpaceType(D);
61 /* check if all function spaces are the same */
62 if (funcspace==FINLEY_ELEMENTS) {
63 reducedIntegrationOrder=FALSE;
64 } else if (funcspace==FINLEY_FACE_ELEMENTS) {
65 reducedIntegrationOrder=FALSE;
66 } else if (funcspace==FINLEY_REDUCED_ELEMENTS) {
67 reducedIntegrationOrder=TRUE;
68 } else if (funcspace==FINLEY_REDUCED_FACE_ELEMENTS) {
69 reducedIntegrationOrder=TRUE;
70 } else {
71 Finley_setError(TYPE_ERROR,"Assemble_LumpedSystem: assemblage failed because of illegal function space.");
72 }
73 if (! Finley_noError()) return;
74
75 /* set all parameters in p*/
76 Assemble_getAssembleParameters(nodes,elements,NULL,lumpedMat, reducedIntegrationOrder, &p);
77 if (! Finley_noError()) return;
78
79 /* check if all function spaces are the same */
80
81 if (! numSamplesEqual(D,p.numQuad,elements->numElements) ) {
82 sprintf(error_msg,"Assemble_LumpedSystem: sample points of coefficient D don't match (%d,%d)",p.numQuad,elements->numElements);
83 Finley_setError(TYPE_ERROR,error_msg);
84 }
85
86 /* check the dimensions: */
87
88 if (p.numEqu==1) {
89 if (!isEmpty(D)) {
90 if (!isDataPointShapeEqual(D,0,dimensions)) {
91 Finley_setError(TYPE_ERROR,"Assemble_LumpedSystem: coefficient D, rank 0 expected.");
92 }
93
94 }
95 } else {
96 if (!isEmpty(D)) {
97 dimensions[0]=p.numEqu;
98 if (!isDataPointShapeEqual(D,1,dimensions)) {
99 sprintf(error_msg,"Assemble_LumpedSystem: coefficient D, expected shape (%d,)",dimensions[0]);
100 Finley_setError(TYPE_ERROR,error_msg);
101 }
102 }
103 }
104
105 if (Finley_noError()) {
106 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 #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp, m_t, diagS)
113 {
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 #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 for (s=0;s<p.row_NS;s++) {
138 rtmp=0;
139 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 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 #endif
156 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 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 } /* end color check */
202 } /* end element loop */
203 } /* end color loop */
204 }
205 } else {
206 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 #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 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 #endif
248 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 /* open loop over all elements: */
255 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 #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 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 #endif
293 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 }
300 } /* end of pointer check */
301 THREAD_MEMFREE(EM_lumpedMat);
302 THREAD_MEMFREE(row_index);
303 } /* end parallel region */
304 }
305 #ifdef Finley_TRACE
306 printf("timing: assemblage lumped PDE: %.4e sec\n",Finley_timer()-time0);
307 #endif
308 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26