/[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 3144 - (show annotations)
Fri Sep 3 00:49:02 2010 UTC (8 years, 7 months ago) by jfenwick
File MIME type: text/plain
File size: 14221 byte(s)
row_node, col_node (Assemble params), node_selection 
removed
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
32 /* Disabled until the tests pass */
33 /* #define NEW_LUMPING */
34
35 /**************************************************************/
36
37 void Dudley_Assemble_LumpedSystem(Dudley_NodeFile* nodes,Dudley_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 #ifdef NEW_LUMPING
51 register double m_t=0., diagS=0.;
52 #endif
53
54 Dudley_resetError();
55
56 if (nodes==NULL || elements==NULL) return;
57 if (isEmpty(lumpedMat) || isEmpty(D)) return;
58 if (isEmpty(lumpedMat) && !isEmpty(D)) {
59 Dudley_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==DUDLEY_ELEMENTS) {
65 reducedIntegrationOrder=FALSE;
66 } else if (funcspace==DUDLEY_FACE_ELEMENTS) {
67 reducedIntegrationOrder=FALSE;
68 } else if (funcspace==DUDLEY_REDUCED_ELEMENTS) {
69 reducedIntegrationOrder=TRUE;
70 } else if (funcspace==DUDLEY_REDUCED_FACE_ELEMENTS) {
71 reducedIntegrationOrder=TRUE;
72 } else {
73 Dudley_setError(TYPE_ERROR,"Assemble_LumpedSystem: assemblage failed because of illegal function space.");
74 }
75 if (! Dudley_noError()) return;
76
77 /* set all parameters in p*/
78 Assemble_getAssembleParameters(nodes,elements,NULL,lumpedMat, reducedIntegrationOrder, &p);
79 if (! Dudley_noError()) return;
80
81 /* check if all function spaces are the same */
82 if (! numSamplesEqual(D,p.numQuadTotal,elements->numElements) ) {
83 sprintf(error_msg,"Assemble_LumpedSystem: sample points of coefficient D don't match (%d,%d)",p.numQuadTotal,elements->numElements);
84 Dudley_setError(TYPE_ERROR,error_msg);
85 }
86
87 /* check the dimensions: */
88 if (p.numEqu==1) {
89 if (!isEmpty(D)) {
90 if (!isDataPointShapeEqual(D,0,dimensions)) {
91 Dudley_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 Dudley_setError(TYPE_ERROR,error_msg);
101 }
102 }
103 }
104 if (Dudley_noError()) {
105 requireWrite(lumpedMat);
106 lumpedMat_p=getSampleDataRW(lumpedMat,0);
107 len_EM_lumpedMat=p.row_numShapesTotal*p.numEqu;
108 len_EM_lumpedMat_size=len_EM_lumpedMat*sizeof(double);
109
110 expandedD=isExpanded(D);
111 S=p.row_jac->BasisFunctions->S;
112
113 #ifdef NEW_LUMPING
114 #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp, diagS, m_t)
115 #else
116 #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp)
117 #endif
118 {
119 EM_lumpedMat=THREAD_MEMALLOC(len_EM_lumpedMat,double);
120 row_index=THREAD_MEMALLOC(p.row_numShapesTotal,index_t);
121 if ( !Dudley_checkPtr(EM_lumpedMat) && !Dudley_checkPtr(row_index) ) {
122 if (p.numEqu == 1) { /* single equation */
123 if (expandedD) { /* with expanded D */
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[INDEX3(0,0,e, p.numQuadTotal,1)]);
130 D_p=getSampleDataRO(D,e);
131 #ifdef NEW_LUMPING /* HRZ lumping */
132 m_t=0; /* mass of the element: m_t */
133 for (q=0;q<p.numQuadTotal;q++) m_t+=Vol[q]*D_p[INDEX2(q, 0,p.numQuadTotal) ];
134
135 diagS=0; /* diagonal sum: S */
136 for (s=0;s<p.row_numShapes;s++) {
137 rtmp=0;
138 for (q=0;q<p.numQuadTotal;q++) rtmp+=Vol[q]*D_p[INDEX2(q, 0, p.numQuadTotal)]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];
139 EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;
140 diagS+=rtmp;
141 }
142 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
143 rtmp=m_t/diagS;
144 for (s=0;s<p.row_numShapes;s++) EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=rtmp;
145
146 #else /* row-sum lumping */
147 for (s=0;s<p.row_numShapes;s++)
148 {
149 rtmp=0;
150 for (q=0;q<p.numQuadTotal;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*D_p[INDEX2(q, 0, p.numQuadTotal)];
151 EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;
152 }
153 #endif
154 for (q=0;q<p.row_numShapesTotal;q++)
155 {
156 row_index[q]=p.row_DOF[elements->Nodes[INDEX2(q,e,p.NN)]];
157 }
158 Dudley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
159 } /* end color check */
160 } /* end element loop */
161 } /* end color loop */
162 } else { /* with constant D */
163
164 for (color=elements->minColor;color<=elements->maxColor;color++) {
165 /* open loop over all elements: */
166 #pragma omp for private(e) schedule(static)
167 for(e=0;e<elements->numElements;e++){
168 if (elements->Color[e]==color) {
169 Vol=&(p.row_jac->volume[INDEX3(0, 0, e, p.numQuadTotal,1)]);
170 D_p=getSampleDataRO(D,e);
171 #ifdef NEW_LUMPING /* HRZ lumping */
172 m_t=0; /* mass of the element: m_t */
173 for (q=0;q<p.numQuadTotal;q++) m_t+=Vol[q];
174 diagS=0; /* diagonal sum: S */
175 for (s=0;s<p.row_numShapes;s++) {
176 rtmp=0;
177 for (q=0;q<p.numQuadTotal;q++)
178 {
179 rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];
180 }
181 EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;
182 diagS+=rtmp;
183 }
184 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
185 rtmp=m_t/diagS*D_p[0];
186 for (s=0;s<p.row_numShapes;s++) EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=rtmp;
187 #else /* row-sum lumping */
188 for (s=0;s<p.row_numShapes;s++) {
189 rtmp=0;
190 for (q=0;q<p.numQuadTotal;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];
191 EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp*D_p[0];
192 }
193 #endif
194 for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(q,e,p.NN)]];
195 Dudley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
196 } /* end color check */
197 } /* end element loop */
198 } /* end color loop */
199
200 }
201 } else { /* system of equation */
202 if (expandedD) { /* with expanded D */
203 for (color=elements->minColor;color<=elements->maxColor;color++) {
204 /* open loop over all elements: */
205 #pragma omp for private(e) schedule(static)
206 for(e=0;e<elements->numElements;e++){
207 if (elements->Color[e]==color) {
208 Vol=&(p.row_jac->volume[INDEX3(0,0,e, p.numQuadTotal,1)]);
209 D_p=getSampleDataRO(D,e);
210
211 #ifdef NEW_LUMPING /* HRZ lumping */
212 for (k=0;k<p.numEqu;k++) {
213 m_t=0; /* mass of the element: m_t */
214 for (q=0;q<p.numQuadTotal;q++) m_t+=Vol[q]*D_p[INDEX3(k,q,0,p.numEqu,p.numQuadTotal)];
215
216 diagS=0; /* diagonal sum: S */
217 for (s=0;s<p.row_numShapes;s++) {
218 rtmp=0;
219 for (q=0;q<p.numQuadTotal;q++) rtmp+=Vol[q]*D_p[INDEX3(k,q,0,p.numEqu,p.numQuadTotal)]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];
220 EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
221 diagS+=rtmp;
222 }
223 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
224 rtmp=m_t/diagS;
225 for (s=0;s<p.row_numShapes;s++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=rtmp;
226 }
227 #else /* row-sum lumping */
228 for (s=0;s<p.row_numShapes;s++) {
229 for (k=0;k<p.numEqu;k++) {
230 rtmp=0.;
231 for (q=0;q<p.numQuadTotal;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*D_p[INDEX3(k,q,0,p.numEqu,p.numQuadTotal)];
232 EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
233 }
234 }
235 #endif
236 for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(q,e,p.NN)]];
237 Dudley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
238 } /* end color check */
239 } /* end element loop */
240 } /* end color loop */
241 } else { /* with constant D */
242 for (color=elements->minColor;color<=elements->maxColor;color++) {
243 /* open loop over all elements: */
244 #pragma omp for private(e) schedule(static)
245 for(e=0;e<elements->numElements;e++){
246 if (elements->Color[e]==color) {
247 Vol=&(p.row_jac->volume[INDEX3(0,0,e, p.numQuadTotal,1)]);
248 D_p=getSampleDataRO(D,e);
249
250 #ifdef NEW_LUMPING /* HRZ lumping */
251 m_t=0; /* mass of the element: m_t */
252 for (q=0;q<p.numQuadTotal;q++) m_t+=Vol[q];
253 diagS=0; /* diagonal sum: S */
254 for (s=0;s<p.row_numShapes;s++) {
255 rtmp=0;
256 for (q=0;q<p.numQuadTotal;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];
257 for (k=0;k<p.numEqu;k++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
258 diagS+=rtmp;
259 }
260 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
261 rtmp=m_t/diagS;
262 for (s=0;s<p.row_numShapes;s++) {
263 for (k=0;k<p.numEqu;k++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=rtmp*D_p[k];
264 }
265 #else /* row-sum lumping */
266 for (s=0;s<p.row_numShapes;s++) {
267 for (k=0;k<p.numEqu;k++) {
268 rtmp=0.;
269 for (q=0;q<p.numQuadTotal;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];
270 EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp*D_p[k];
271 }
272 }
273 #endif
274 for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(q,e,p.NN)]];
275 Dudley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
276 } /* end color check */
277 } /* end element loop */
278 } /* end color loop */
279 }
280 }
281 } /* end of pointer check */
282 THREAD_MEMFREE(EM_lumpedMat);
283 THREAD_MEMFREE(row_index);
284 } /* end parallel region */
285 }
286 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26