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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3136 - (hide annotations)
Thu Sep 2 03:48:29 2010 UTC (8 years, 11 months ago) by jfenwick
File MIME type: text/plain
File size: 14350 byte(s)
Removed numSubElements
1 jgs 82
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4 jfenwick 2881 * Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1811 * 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 ksteube 1312
14 ksteube 1811
15 jgs 82 /**************************************************************/
16    
17 gross 1204 /* assembles the mass matrix in lumped form */
18 jgs 82
19 gross 1204 /* The coefficient D has to be defined on the integration points or not present. */
20 jgs 82
21 gross 1204 /* lumpedMat has to be initialized before the routine is called. */
22 jgs 82
23     /**************************************************************/
24    
25     #include "Assemble.h"
26     #include "Util.h"
27     #ifdef _OPENMP
28     #include <omp.h>
29     #endif
30    
31    
32 ksteube 1221 /* Disabled until the tests pass */
33 gross 2989 /* #define NEW_LUMPING */
34 btully 1220
35 jgs 82 /**************************************************************/
36    
37 jfenwick 3086 void Dudley_Assemble_LumpedSystem(Dudley_NodeFile* nodes,Dudley_ElementFile* elements, escriptDataC* lumpedMat, escriptDataC* D)
38 gross 1204 {
39 jgs 82
40 gross 1204 bool_t reducedIntegrationOrder=FALSE, expandedD;
41 jgs 150 char error_msg[LenErrorMsg_MAX];
42 gross 798 Assemble_Parameters p;
43 jfenwick 3136 dim_t dimensions[ESCRIPT_MAX_DATA_RANK], k, e, len_EM_lumpedMat, q, s;
44 gross 1028 type_t funcspace;
45 gross 1204 index_t color,*row_index=NULL;
46 jfenwick 2271 __const double *D_p=NULL;
47     double *S=NULL, *EM_lumpedMat=NULL, *Vol=NULL, *lumpedMat_p=NULL;
48 phornby 1628 register double rtmp;
49 gross 1204 size_t len_EM_lumpedMat_size;
50 caltinay 2750 #ifdef NEW_LUMPING
51 gross 2748 register double m_t=0., diagS=0.;
52 caltinay 2750 #endif
53 gross 2989
54 jfenwick 3086 Dudley_resetError();
55 jgs 82
56     if (nodes==NULL || elements==NULL) return;
57 gross 1204 if (isEmpty(lumpedMat) || isEmpty(D)) return;
58     if (isEmpty(lumpedMat) && !isEmpty(D)) {
59 jfenwick 3086 Dudley_setError(TYPE_ERROR,"Assemble_LumpedSystem: coefficients are non-zero but no lumped matrix is given.");
60 gross 1204 return;
61 gross 798 }
62 gross 1204 funcspace=getFunctionSpaceType(D);
63 jgs 82 /* check if all function spaces are the same */
64 jfenwick 3086 if (funcspace==DUDLEY_ELEMENTS) {
65 gross 798 reducedIntegrationOrder=FALSE;
66 jfenwick 3086 } else if (funcspace==DUDLEY_FACE_ELEMENTS) {
67 gross 798 reducedIntegrationOrder=FALSE;
68 jfenwick 3086 } else if (funcspace==DUDLEY_REDUCED_ELEMENTS) {
69 gross 1072 reducedIntegrationOrder=TRUE;
70 jfenwick 3086 } else if (funcspace==DUDLEY_REDUCED_FACE_ELEMENTS) {
71 gross 1072 reducedIntegrationOrder=TRUE;
72 gross 798 } else {
73 jfenwick 3086 Dudley_setError(TYPE_ERROR,"Assemble_LumpedSystem: assemblage failed because of illegal function space.");
74 gross 798 }
75 jfenwick 3086 if (! Dudley_noError()) return;
76 jgs 82
77 gross 798 /* set all parameters in p*/
78 gross 1204 Assemble_getAssembleParameters(nodes,elements,NULL,lumpedMat, reducedIntegrationOrder, &p);
79 jfenwick 3086 if (! Dudley_noError()) return;
80 gross 2989
81 gross 798 /* check if all function spaces are the same */
82 gross 2989 if (! numSamplesEqual(D,p.numQuadTotal,elements->numElements) ) {
83 gross 2748 sprintf(error_msg,"Assemble_LumpedSystem: sample points of coefficient D don't match (%d,%d)",p.numQuadSub,elements->numElements);
84 jfenwick 3086 Dudley_setError(TYPE_ERROR,error_msg);
85 jgs 82 }
86    
87     /* check the dimensions: */
88 gross 1204 if (p.numEqu==1) {
89 jgs 82 if (!isEmpty(D)) {
90     if (!isDataPointShapeEqual(D,0,dimensions)) {
91 jfenwick 3086 Dudley_setError(TYPE_ERROR,"Assemble_LumpedSystem: coefficient D, rank 0 expected.");
92 jgs 82 }
93 gross 1204
94 jgs 82 }
95     } else {
96     if (!isEmpty(D)) {
97     dimensions[0]=p.numEqu;
98 gross 1204 if (!isDataPointShapeEqual(D,1,dimensions)) {
99     sprintf(error_msg,"Assemble_LumpedSystem: coefficient D, expected shape (%d,)",dimensions[0]);
100 jfenwick 3086 Dudley_setError(TYPE_ERROR,error_msg);
101 jgs 82 }
102     }
103     }
104 jfenwick 3086 if (Dudley_noError()) {
105 jfenwick 2271 requireWrite(lumpedMat);
106     lumpedMat_p=getSampleDataRW(lumpedMat,0);
107 gross 2748 len_EM_lumpedMat=p.row_numShapesTotal*p.numEqu;
108 gross 1204 len_EM_lumpedMat_size=len_EM_lumpedMat*sizeof(double);
109 gross 2748
110 gross 1204 expandedD=isExpanded(D);
111 gross 2748 S=p.row_jac->BasisFunctions->S;
112 gross 2989
113 caltinay 2750 #ifdef NEW_LUMPING
114 jfenwick 3136 #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp, diagS, m_t)
115 caltinay 2750 #else
116 jfenwick 3136 #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp)
117 caltinay 2750 #endif
118 gross 1204 {
119     EM_lumpedMat=THREAD_MEMALLOC(len_EM_lumpedMat,double);
120 gross 2748 row_index=THREAD_MEMALLOC(p.row_numShapesTotal,index_t);
121 jfenwick 3086 if ( !Dudley_checkPtr(EM_lumpedMat) && !Dudley_checkPtr(row_index) ) {
122 gross 2988 if (p.numEqu == 1) { /* single equation */
123     if (expandedD) { /* with expanded D */
124 gross 1204 for (color=elements->minColor;color<=elements->maxColor;color++) {
125     /* open loop over all elements: */
126     #pragma omp for private(e) schedule(static)
127 gross 2988 for(e=0;e<elements->numElements;e++){
128     if (elements->Color[e]==color) {
129 jfenwick 3136 Vol=&(p.row_jac->volume[INDEX3(0,0,e, p.numQuadSub,1)]);
130 jfenwick 2770 D_p=getSampleDataRO(D,e);
131 gross 2988 #ifdef NEW_LUMPING /* HRZ lumping */
132 gross 2748 m_t=0; /* mass of the element: m_t */
133 jfenwick 3136 for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q]*D_p[INDEX2(q, 0,p.numQuadSub) ];
134 gross 2988
135 gross 2748 diagS=0; /* diagonal sum: S */
136     for (s=0;s<p.row_numShapes;s++) {
137     rtmp=0;
138 jfenwick 3136 for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*D_p[INDEX2(q, 0, p.numQuadSub)]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];
139 gross 2748 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 gross 2988 rtmp=m_t/diagS;
144 gross 2748 for (s=0;s<p.row_numShapes;s++) EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=rtmp;
145    
146     #else /* row-sum lumping */
147 jfenwick 3136 for (s=0;s<p.row_numShapes;s++)
148     {
149 gross 2748 rtmp=0;
150 jfenwick 3136 for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*D_p[INDEX2(q, 0, p.numQuadSub)];
151 gross 2989 EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;
152 gross 2748 }
153     #endif
154 jfenwick 3136 for (q=0;q<p.row_numShapesTotal;q++)
155     {
156     row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[INDEX2(q,0, p.row_numShapesTotal)],e,p.NN)]];
157     }
158 jfenwick 3086 Dudley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
159 gross 2988 } /* end color check */
160 gross 1204 } /* end element loop */
161     } /* end color loop */
162 gross 2988 } else { /* with constant D */
163 gross 2990
164 gross 1204 for (color=elements->minColor;color<=elements->maxColor;color++) {
165 gross 2748 /* open loop over all elements: */
166     #pragma omp for private(e) schedule(static)
167     for(e=0;e<elements->numElements;e++){
168 gross 2988 if (elements->Color[e]==color) {
169 jfenwick 3136 Vol=&(p.row_jac->volume[INDEX3(0, 0, e, p.numQuadSub,1)]);
170 jfenwick 2770 D_p=getSampleDataRO(D,e);
171 gross 2988 #ifdef NEW_LUMPING /* HRZ lumping */
172 gross 2748 m_t=0; /* mass of the element: m_t */
173     for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q];
174     diagS=0; /* diagonal sum: S */
175     for (s=0;s<p.row_numShapes;s++) {
176     rtmp=0;
177 jfenwick 3136 for (q=0;q<p.numQuadSub;q++)
178     {
179     rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];
180     }
181 gross 2748 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 gross 2988 rtmp=m_t/diagS*D_p[0];
186 gross 2748 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.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];
191 gross 2988 EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp*D_p[0];
192 gross 2748 }
193     #endif
194 jfenwick 3136 for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[INDEX2(q, 0, p.row_numShapesTotal)],e,p.NN)]];
195 jfenwick 3086 Dudley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
196 gross 2748 } /* end color check */
197 gross 1204 } /* end element loop */
198 gross 2748 } /* end color loop */
199    
200 gross 1204 }
201 gross 2988 } else { /* system of equation */
202     if (expandedD) { /* with expanded D */
203 gross 1204 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 jfenwick 3136 Vol=&(p.row_jac->volume[INDEX3(0,0,e, p.numQuadSub,1)]);
209 gross 2988 D_p=getSampleDataRO(D,e);
210    
211 gross 2748 #ifdef NEW_LUMPING /* HRZ lumping */
212     for (k=0;k<p.numEqu;k++) {
213     m_t=0; /* mass of the element: m_t */
214 jfenwick 3136 for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q]*D_p[INDEX3(k,q,0,p.numEqu,p.numQuadSub)];
215 gross 2748
216     diagS=0; /* diagonal sum: S */
217     for (s=0;s<p.row_numShapes;s++) {
218     rtmp=0;
219 jfenwick 3136 for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*D_p[INDEX3(k,q,0,p.numEqu,p.numQuadSub)]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];
220 gross 2748 EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
221     diagS+=rtmp;
222 gross 2988 }
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 jfenwick 3136 for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*D_p[INDEX3(k,q,0,p.numEqu,p.numQuadSub)];
232 gross 2988 EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
233     }
234     }
235     #endif
236 jfenwick 3136 for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[INDEX2(q,0,p.row_numShapesTotal)],e,p.NN)]];
237 jfenwick 3086 Dudley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
238 gross 1204 } /* end color check */
239     } /* end element loop */
240     } /* end color loop */
241 gross 2988 } else { /* with constant D */
242 gross 1204 for (color=elements->minColor;color<=elements->maxColor;color++) {
243 gross 2748 /* open loop over all elements: */
244 gross 1204 #pragma omp for private(e) schedule(static)
245     for(e=0;e<elements->numElements;e++){
246     if (elements->Color[e]==color) {
247 jfenwick 3136 Vol=&(p.row_jac->volume[INDEX3(0,0,e, p.numQuadSub,1)]);
248 jfenwick 2770 D_p=getSampleDataRO(D,e);
249 gross 2748
250     #ifdef NEW_LUMPING /* HRZ lumping */
251     m_t=0; /* mass of the element: m_t */
252 gross 2988 for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q];
253 gross 2748 diagS=0; /* diagonal sum: S */
254     for (s=0;s<p.row_numShapes;s++) {
255     rtmp=0;
256     for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];
257 gross 2988 for (k=0;k<p.numEqu;k++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
258 gross 2748 diagS+=rtmp;
259     }
260 gross 2988 /* 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 gross 2748 for (s=0;s<p.row_numShapes;s++) {
267 gross 2988 for (k=0;k<p.numEqu;k++) {
268 gross 2748 rtmp=0.;
269     for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];
270 gross 2989 EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp*D_p[k];
271 gross 2748 }
272 gross 2988 }
273     #endif
274 jfenwick 3136 for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[INDEX2(q,0,p.row_numShapesTotal)],e,p.NN)]];
275 jfenwick 3086 Dudley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
276 gross 1204 } /* end color check */
277     } /* end element loop */
278     } /* end color loop */
279     }
280 gross 798 }
281 gross 1204 } /* end of pointer check */
282     THREAD_MEMFREE(EM_lumpedMat);
283     THREAD_MEMFREE(row_index);
284     } /* end parallel region */
285 jgs 82 }
286     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26