/[escript]/branches/windows_from_1456_trunk_1574_merged_in/finley/src/Assemble_LumpedSystem.c
ViewVC logotype

Contents of /branches/windows_from_1456_trunk_1574_merged_in/finley/src/Assemble_LumpedSystem.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1580 - (show annotations)
Tue Jun 3 14:03:40 2008 UTC (13 years, 5 months ago) by phornby
File MIME type: text/plain
File size: 14882 byte(s)
This once again compiles and links under windows after merging to trunk version 1574. scons run_tests also passes now. py_tests still fails in the same way.  Had to link against winsock2 to get gethostname(), and had to re-do some work on eliminating unused vars. Also eliminated signed/unsigned comparisons where I saw them.
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 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 double *S=NULL, *EM_lumpedMat=NULL, *Vol=NULL, *D_p=NULL, *lumpedMat_p=NULL;
48 register double rtmp;
49 size_t len_EM_lumpedMat_size;
50 #ifdef NEW_LUMPING /* HRZ 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, m_t, diagS)
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 m_t=0; /* mass of the element: m_t */
135 for (q=0;q<p.numQuad;q++) {
136 m_t+=Vol[q]*D_p[q];
137 }
138 diagS=0; /* diagonal sum: S */
139 for (s=0;s<p.row_NS;s++) {
140 rtmp=0;
141 for (q=0;q<p.numQuad;q++) {
142 rtmp+=Vol[q]*D_p[q]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
143 }
144 EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
145 diagS+=EM_lumpedMat[INDEX2(0,s,p.numEqu)];
146 }
147 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
148 for (s=0;s<p.row_NS;s++) {
149 EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=m_t/diagS;
150 }
151 #else /* row-sum lumping */
152 for (s=0;s<p.row_NS;s++) {
153 rtmp=0;
154 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[q];
155 EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
156 }
157 #endif
158 for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
159 Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
160 } /* end color check */
161 } /* end element loop */
162 } /* end color loop */
163 } else {
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[INDEX2(0,e,p.numQuad)]);
170 memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
171 D_p=getSampleData(D,e);
172 #ifdef NEW_LUMPING /* HRZ lumping */
173 /*
174 * Number of PDEs: 1
175 * D_p varies over element: False
176 */
177 m_t=0; /* mass of the element: m_t */
178 for (q=0;q<p.numQuad;q++) {
179 m_t+=Vol[q]*D_p[0];
180 }
181 diagS=0; /* diagonal sum: S */
182 for (s=0;s<p.row_NS;s++) {
183 rtmp=0;
184 for (q=0;q<p.numQuad;q++) {
185 rtmp+=Vol[q]*D_p[0]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
186 }
187 EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
188 diagS+=EM_lumpedMat[INDEX2(0,s,p.numEqu)];
189 }
190 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
191 for (s=0;s<p.row_NS;s++) {
192 EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=m_t/diagS;
193 }
194 #else /* row-sum lumping */
195 for (s=0;s<p.row_NS;s++) {
196 rtmp=0;
197 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
198 EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp*D_p[0];
199 }
200 #endif
201 for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
202 Finley_Util_AddScatter(p.row_NN,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 } else {
208 if (expandedD) {
209 for (color=elements->minColor;color<=elements->maxColor;color++) {
210 /* open loop over all elements: */
211 #pragma omp for private(e) schedule(static)
212 for(e=0;e<elements->numElements;e++){
213 if (elements->Color[e]==color) {
214 Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
215 memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
216 D_p=getSampleData(D,e);
217 #ifdef NEW_LUMPING /* HRZ lumping */
218 /*
219 * Number of PDEs: Multiple
220 * D_p varies over element: True
221 */
222 for (k=0;k<p.numEqu;k++) {
223 m_t=0; /* mass of the element: m_t */
224 for (q=0;q<p.numQuad;q++) {
225 m_t+=Vol[q]*D_p[INDEX2(k,q,p.numEqu)];
226 }
227 diagS=0; /* diagonal sum: S */
228 for (s=0;s<p.row_NS;s++) {
229 rtmp=0;
230 for (q=0;q<p.numQuad;q++) {
231 rtmp+=Vol[q]*D_p[INDEX2(k,q,p.numEqu)]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
232 }
233 EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
234 diagS+=EM_lumpedMat[INDEX2(k,s,p.numEqu)];
235 }
236 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
237 for (s=0;s<p.row_NS;s++) {
238 EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=m_t/diagS;
239 }
240 }
241 #else /* row-sum lumping */
242 for (s=0;s<p.row_NS;s++) {
243 for (k=0;k<p.numEqu;k++) {
244 rtmp=0.;
245 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[INDEX2(k,q,p.numEqu)];
246 EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
247 }
248 }
249 #endif
250 for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
251 Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
252 } /* end color check */
253 } /* end element loop */
254 } /* end color loop */
255 } else {
256 /* open loop over all elements: */
257 for (color=elements->minColor;color<=elements->maxColor;color++) {
258 #pragma omp for private(e) schedule(static)
259 for(e=0;e<elements->numElements;e++){
260 if (elements->Color[e]==color) {
261 Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
262 memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
263 D_p=getSampleData(D,e);
264 #ifdef NEW_LUMPING /* HRZ lumping */
265 /*
266 * Number of PDEs: Multiple
267 * D_p varies over element: False
268 */
269 for (k=0;k<p.numEqu;k++) {
270 m_t=0; /* mass of the element: m_t */
271 for (q=0;q<p.numQuad;q++) {
272 m_t+=Vol[q]*D_p[k];
273 }
274 diagS=0; /* diagonal sum: S */
275 for (s=0;s<p.row_NS;s++) {
276 rtmp=0;
277 for (q=0;q<p.numQuad;q++) {
278 rtmp+=Vol[q]*D_p[k]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
279 }
280 EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
281 diagS+=EM_lumpedMat[INDEX2(k,s,p.numEqu)];
282 }
283 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
284 for (s=0;s<p.row_NS;s++) {
285 EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=m_t/diagS;
286 }
287 }
288 #else /* row-sum lumping */
289 for (s=0;s<p.row_NS;s++) {
290 rtmp=0;
291 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
292 for (k=0;k<p.numEqu;k++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp*D_p[k];
293 }
294 #endif
295 for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
296 Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
297 } /* end color check */
298 } /* end element loop */
299 } /* end color loop */
300 }
301 }
302 } /* end of pointer check */
303 THREAD_MEMFREE(EM_lumpedMat);
304 THREAD_MEMFREE(row_index);
305 } /* end parallel region */
306 }
307 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26