/[escript]/branches/domexper/dudley/src/Assemble_LumpedSystem.c
ViewVC logotype

Contents of /branches/domexper/dudley/src/Assemble_LumpedSystem.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3205 - (show annotations)
Fri Sep 24 00:30:43 2010 UTC (8 years, 8 months ago) by jfenwick
File MIME type: text/plain
File size: 14048 byte(s)
Removing -total from field names

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 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 #include "ShapeTable.h"
32
33 /* Disabled until the tests pass */
34 /* #define NEW_LUMPING */
35
36 /**************************************************************/
37
38 void Dudley_Assemble_LumpedSystem(Dudley_NodeFile* nodes,Dudley_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 __const double *D_p=NULL;
48 const double* S=NULL;
49 double *EM_lumpedMat=NULL, *lumpedMat_p=NULL;
50 register double rtmp;
51 size_t len_EM_lumpedMat_size;
52 #ifdef NEW_LUMPING
53 register double m_t=0., diagS=0.;
54 #endif
55
56 Dudley_resetError();
57
58 if (nodes==NULL || elements==NULL) return;
59 if (isEmpty(lumpedMat) || isEmpty(D)) return;
60 if (isEmpty(lumpedMat) && !isEmpty(D)) {
61 Dudley_setError(TYPE_ERROR,"Assemble_LumpedSystem: coefficients are non-zero but no lumped matrix is given.");
62 return;
63 }
64 funcspace=getFunctionSpaceType(D);
65 /* check if all function spaces are the same */
66 if (funcspace==DUDLEY_ELEMENTS) {
67 reducedIntegrationOrder=FALSE;
68 } else if (funcspace==DUDLEY_FACE_ELEMENTS) {
69 reducedIntegrationOrder=FALSE;
70 } else if (funcspace==DUDLEY_REDUCED_ELEMENTS) {
71 reducedIntegrationOrder=TRUE;
72 } else if (funcspace==DUDLEY_REDUCED_FACE_ELEMENTS) {
73 reducedIntegrationOrder=TRUE;
74 } else {
75 Dudley_setError(TYPE_ERROR,"Assemble_LumpedSystem: assemblage failed because of illegal function space.");
76 }
77 if (! Dudley_noError()) return;
78
79 /* set all parameters in p*/
80 Assemble_getAssembleParameters(nodes,elements,NULL,lumpedMat, reducedIntegrationOrder, &p);
81 if (! Dudley_noError()) return;
82
83 /* check if all function spaces are the same */
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 Dudley_setError(TYPE_ERROR,error_msg);
87 }
88
89 /* check the dimensions: */
90 if (p.numEqu==1) {
91 if (!isEmpty(D)) {
92 if (!isDataPointShapeEqual(D,0,dimensions)) {
93 Dudley_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 Dudley_setError(TYPE_ERROR,error_msg);
103 }
104 }
105 }
106 if (Dudley_noError()) {
107 requireWrite(lumpedMat);
108 lumpedMat_p=getSampleDataRW(lumpedMat,0);
109 len_EM_lumpedMat=p.numShapes*p.numEqu;
110 len_EM_lumpedMat_size=len_EM_lumpedMat*sizeof(double);
111
112 expandedD=isExpanded(D);
113 // S=p.row_jac->BasisFunctions->S;
114 if (!getQuadShape(elements->numDim, reducedIntegrationOrder, &S))
115 {
116 Dudley_setError(TYPE_ERROR, "Assemble_LumpedSystem: Unable to locate shape function.");
117 }
118
119 #ifdef NEW_LUMPING
120 #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp, diagS, m_t)
121 #else
122 #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp)
123 #endif
124 {
125 EM_lumpedMat=THREAD_MEMALLOC(len_EM_lumpedMat,double);
126 row_index=THREAD_MEMALLOC(p.numShapes,index_t);
127 if ( !Dudley_checkPtr(EM_lumpedMat) && !Dudley_checkPtr(row_index) ) {
128 if (p.numEqu == 1) { /* single equation */
129 if (expandedD) { /* with expanded D */
130 for (color=elements->minColor;color<=elements->maxColor;color++) {
131 /* open loop over all elements: */
132 #pragma omp for private(e) schedule(static)
133 for(e=0;e<elements->numElements;e++){
134 if (elements->Color[e]==color) {
135 // Vol=&(p.row_jac->volume[INDEX3(0,0,e, p.numQuad,1)]);
136 double vol=p.row_jac->absD[e]*p.row_jac->quadweight;
137 D_p=getSampleDataRO(D,e);
138 #ifdef NEW_LUMPING /* HRZ lumping */
139 m_t=0; /* mass of the element: m_t */
140 for (q=0;q<p.numQuad;q++) m_t+=vol*D_p[INDEX2(q, 0,p.numQuad) ];
141
142 diagS=0; /* diagonal sum: S */
143 for (s=0;s<p.numShapes;s++) {
144 rtmp=0;
145 for (q=0;q<p.numQuad;q++) rtmp+=vol*D_p[INDEX2(q, 0, p.numQuad)]*S[INDEX2(s,q,p.numShapes)]*S[INDEX2(s,q,p.numShapes)];
146 EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;
147 diagS+=rtmp;
148 }
149 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
150 rtmp=m_t/diagS;
151 for (s=0;s<p.numShapes;s++) EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=rtmp;
152
153 #else /* row-sum lumping */
154 for (s=0;s<p.numShapes;s++)
155 {
156 rtmp=0;
157 for (q=0;q<p.numQuad;q++) rtmp+=vol*S[INDEX2(s,q,p.numShapes)]*D_p[INDEX2(q, 0, p.numQuad)];
158 EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;
159 }
160 #endif
161 for (q=0;q<p.numShapes;q++)
162 {
163 row_index[q]=p.row_DOF[elements->Nodes[INDEX2(q,e,p.NN)]];
164 }
165 Dudley_Util_AddScatter(p.numShapes,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 { /* with constant D */
170
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 double vol=p.row_jac->absD[e]*p.row_jac->quadweight;
177 D_p=getSampleDataRO(D,e);
178 #ifdef NEW_LUMPING /* HRZ lumping */
179 m_t=0; /* mass of the element: m_t */
180 for (q=0;q<p.numQuad;q++) m_t+=vol;
181 diagS=0; /* diagonal sum: S */
182 for (s=0;s<p.numShapes;s++) {
183 rtmp=0;
184 for (q=0;q<p.numQuad;q++)
185 {
186 rtmp+=vol*S[INDEX2(s,q,p.numShapes)]*S[INDEX2(s,q,p.numShapes)];
187 }
188 EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;
189 diagS+=rtmp;
190 }
191 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
192 rtmp=m_t/diagS*D_p[0];
193 for (s=0;s<p.numShapes;s++) EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=rtmp;
194 #else /* row-sum lumping */
195 for (s=0;s<p.numShapes;s++) {
196 rtmp=0;
197 for (q=0;q<p.numQuad;q++) rtmp+=vol*S[INDEX2(s,q,p.numShapes)];
198 EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp*D_p[0];
199 }
200 #endif
201 for (q=0;q<p.numShapes;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(q,e,p.NN)]];
202 Dudley_Util_AddScatter(p.numShapes,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
203 } /* end color check */
204 } /* end element loop */
205 } /* end color loop */
206
207 }
208 } else { /* system of equation */
209 if (expandedD) { /* with expanded D */
210 for (color=elements->minColor;color<=elements->maxColor;color++) {
211 /* open loop over all elements: */
212 #pragma omp for private(e) schedule(static)
213 for(e=0;e<elements->numElements;e++){
214 if (elements->Color[e]==color) {
215 double vol=p.row_jac->absD[e]*p.row_jac->quadweight;
216 D_p=getSampleDataRO(D,e);
217
218 #ifdef NEW_LUMPING /* HRZ lumping */
219 for (k=0;k<p.numEqu;k++) {
220 m_t=0; /* mass of the element: m_t */
221 for (q=0;q<p.numQuad;q++) m_t+=vol*D_p[INDEX3(k,q,0,p.numEqu,p.numQuad)];
222
223 diagS=0; /* diagonal sum: S */
224 for (s=0;s<p.numShapes;s++) {
225 rtmp=0;
226 for (q=0;q<p.numQuad;q++) rtmp+=vol*D_p[INDEX3(k,q,0,p.numEqu,p.numQuad)]*S[INDEX2(s,q,p.numShapes)]*S[INDEX2(s,q,p.numShapes)];
227 EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
228 diagS+=rtmp;
229 }
230 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
231 rtmp=m_t/diagS;
232 for (s=0;s<p.numShapes;s++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=rtmp;
233 }
234 #else /* row-sum lumping */
235 for (s=0;s<p.numShapes;s++) {
236 for (k=0;k<p.numEqu;k++) {
237 rtmp=0.;
238 for (q=0;q<p.numQuad;q++) rtmp+=vol*S[INDEX2(s,q,p.numShapes)]*D_p[INDEX3(k,q,0,p.numEqu,p.numQuad)];
239 EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
240 }
241 }
242 #endif
243 for (q=0;q<p.numShapes;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(q,e,p.NN)]];
244 Dudley_Util_AddScatter(p.numShapes,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
245 } /* end color check */
246 } /* end element loop */
247 } /* end color loop */
248 } else { /* with constant D */
249 for (color=elements->minColor;color<=elements->maxColor;color++) {
250 /* open loop over all elements: */
251 #pragma omp for private(e) schedule(static)
252 for(e=0;e<elements->numElements;e++){
253 if (elements->Color[e]==color) {
254 double vol=p.row_jac->absD[e]*p.row_jac->quadweight; D_p=getSampleDataRO(D,e);
255
256 #ifdef NEW_LUMPING /* HRZ lumping */
257 m_t=0; /* mass of the element: m_t */
258 for (q=0;q<p.numQuad;q++) m_t+=vol;
259 diagS=0; /* diagonal sum: S */
260 for (s=0;s<p.numShapes;s++) {
261 rtmp=0;
262 for (q=0;q<p.numQuad;q++) rtmp+=vol*S[INDEX2(s,q,p.numShapes)]*S[INDEX2(s,q,p.numShapes)];
263 for (k=0;k<p.numEqu;k++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
264 diagS+=rtmp;
265 }
266 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
267 rtmp=m_t/diagS;
268 for (s=0;s<p.numShapes;s++) {
269 for (k=0;k<p.numEqu;k++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=rtmp*D_p[k];
270 }
271 #else /* row-sum lumping */
272 for (s=0;s<p.numShapes;s++) {
273 for (k=0;k<p.numEqu;k++) {
274 rtmp=0.;
275 for (q=0;q<p.numQuad;q++) rtmp+=vol*S[INDEX2(s,q,p.numShapes)];
276 EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp*D_p[k];
277 }
278 }
279 #endif
280 for (q=0;q<p.numShapes;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(q,e,p.NN)]];
281 Dudley_Util_AddScatter(p.numShapes,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
282 } /* end color check */
283 } /* end element loop */
284 } /* end color loop */
285 }
286 }
287 } /* end of pointer check */
288 THREAD_MEMFREE(EM_lumpedMat);
289 THREAD_MEMFREE(row_index);
290 } /* end parallel region */
291 }
292 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26