/[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 1811 - (show annotations)
Thu Sep 25 23:11:13 2008 UTC (10 years, 6 months ago) by ksteube
File MIME type: text/plain
File size: 14940 byte(s)
Copyright updated in all files

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 double *S=NULL, *EM_lumpedMat=NULL, *Vol=NULL, *D_p=NULL, *lumpedMat_p=NULL;
47 register double rtmp;
48 size_t len_EM_lumpedMat_size;
49
50 #ifdef NEW_LUMPING
51 register double m_t, diagS;
52 #endif
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)
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 for (color=elements->minColor;color<=elements->maxColor;color++) {
122 /* open loop over all elements: */
123 #pragma omp for private(e) schedule(static)
124 for(e=0;e<elements->numElements;e++){
125 if (elements->Color[e]==color) {
126 Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
127 memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
128 D_p=getSampleData(D,e);
129 #ifdef NEW_LUMPING /* HRZ lumping */
130 /*
131 * Number of PDEs: 1
132 * D_p varies over element: True
133 */
134 #pragma omp parallel private(m_t)
135 m_t=0; /* mass of the element: m_t */
136 for (q=0;q<p.numQuad;q++) {
137 m_t+=Vol[q]*D_p[q];
138 }
139 #pragma omp parallel private(diagS)
140 diagS=0; /* diagonal sum: S */
141 for (s=0;s<p.row_NS;s++) {
142 rtmp=0;
143 for (q=0;q<p.numQuad;q++) {
144 rtmp+=Vol[q]*D_p[q]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
145 }
146 EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
147 diagS+=EM_lumpedMat[INDEX2(0,s,p.numEqu)];
148 }
149 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
150 for (s=0;s<p.row_NS;s++) {
151 EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=m_t/diagS;
152 }
153 #else /* row-sum lumping */
154 for (s=0;s<p.row_NS;s++) {
155 rtmp=0;
156 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[q];
157 EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
158 }
159 #endif
160 for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
161 Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
162 } /* end color check */
163 } /* end element loop */
164 } /* end color loop */
165 } else {
166 for (color=elements->minColor;color<=elements->maxColor;color++) {
167 /* open loop over all elements: */
168 #pragma omp for private(e) schedule(static)
169 for(e=0;e<elements->numElements;e++){
170 if (elements->Color[e]==color) {
171 Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
172 memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
173 D_p=getSampleData(D,e);
174 #ifdef NEW_LUMPING /* HRZ lumping */
175 /*
176 * Number of PDEs: 1
177 * D_p varies over element: False
178 */
179 m_t=0; /* mass of the element: m_t */
180 for (q=0;q<p.numQuad;q++) {
181 m_t+=Vol[q]*D_p[0];
182 }
183 diagS=0; /* diagonal sum: S */
184 for (s=0;s<p.row_NS;s++) {
185 rtmp=0;
186 for (q=0;q<p.numQuad;q++) {
187 rtmp+=Vol[q]*D_p[0]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
188 }
189 EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
190 diagS+=EM_lumpedMat[INDEX2(0,s,p.numEqu)];
191 }
192 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
193 for (s=0;s<p.row_NS;s++) {
194 EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=m_t/diagS;
195 }
196 #else /* row-sum lumping */
197 for (s=0;s<p.row_NS;s++) {
198 rtmp=0;
199 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
200 EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp*D_p[0];
201 }
202 #endif
203 for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
204 Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
205 } /* end color check */
206 } /* end element loop */
207 } /* end color loop */
208 }
209 } else {
210 if (expandedD) {
211 for (color=elements->minColor;color<=elements->maxColor;color++) {
212 /* open loop over all elements: */
213 #pragma omp for private(e) schedule(static)
214 for(e=0;e<elements->numElements;e++){
215 if (elements->Color[e]==color) {
216 Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
217 memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
218 D_p=getSampleData(D,e);
219 #ifdef NEW_LUMPING /* HRZ lumping */
220 /*
221 * Number of PDEs: Multiple
222 * D_p varies over element: True
223 */
224 for (k=0;k<p.numEqu;k++) {
225 m_t=0; /* mass of the element: m_t */
226 for (q=0;q<p.numQuad;q++) {
227 m_t+=Vol[q]*D_p[INDEX2(k,q,p.numEqu)];
228 }
229 diagS=0; /* diagonal sum: S */
230 for (s=0;s<p.row_NS;s++) {
231 rtmp=0;
232 for (q=0;q<p.numQuad;q++) {
233 rtmp+=Vol[q]*D_p[INDEX2(k,q,p.numEqu)]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
234 }
235 EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
236 diagS+=EM_lumpedMat[INDEX2(k,s,p.numEqu)];
237 }
238 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
239 for (s=0;s<p.row_NS;s++) {
240 EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=m_t/diagS;
241 }
242 }
243 #else /* row-sum lumping */
244 for (s=0;s<p.row_NS;s++) {
245 for (k=0;k<p.numEqu;k++) {
246 rtmp=0.;
247 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[INDEX2(k,q,p.numEqu)];
248 EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
249 }
250 }
251 #endif
252 for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
253 Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
254 } /* end color check */
255 } /* end element loop */
256 } /* end color loop */
257 } else {
258 /* open loop over all elements: */
259 for (color=elements->minColor;color<=elements->maxColor;color++) {
260 #pragma omp for private(e) schedule(static)
261 for(e=0;e<elements->numElements;e++){
262 if (elements->Color[e]==color) {
263 Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
264 memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
265 D_p=getSampleData(D,e);
266 #ifdef NEW_LUMPING /* HRZ lumping */
267 /*
268 * Number of PDEs: Multiple
269 * D_p varies over element: False
270 */
271 for (k=0;k<p.numEqu;k++) {
272 m_t=0; /* mass of the element: m_t */
273 for (q=0;q<p.numQuad;q++) {
274 m_t+=Vol[q]*D_p[k];
275 }
276 diagS=0; /* diagonal sum: S */
277 for (s=0;s<p.row_NS;s++) {
278 rtmp=0;
279 for (q=0;q<p.numQuad;q++) {
280 rtmp+=Vol[q]*D_p[k]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
281 }
282 EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
283 diagS+=EM_lumpedMat[INDEX2(k,s,p.numEqu)];
284 }
285 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
286 for (s=0;s<p.row_NS;s++) {
287 EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=m_t/diagS;
288 }
289 }
290 #else /* row-sum lumping */
291 for (s=0;s<p.row_NS;s++) {
292 rtmp=0;
293 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
294 for (k=0;k<p.numEqu;k++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp*D_p[k];
295 }
296 #endif
297 for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
298 Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
299 } /* end color check */
300 } /* end element loop */
301 } /* end color loop */
302 }
303 }
304 } /* end of pointer check */
305 THREAD_MEMFREE(EM_lumpedMat);
306 THREAD_MEMFREE(row_index);
307 } /* end parallel region */
308 }
309 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26