/[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 1204 - (show annotations)
Sat Jun 23 11:43:12 2007 UTC (11 years, 9 months ago) by gross
File MIME type: text/plain
File size: 10214 byte(s)
a frame for an improved lumping procedure added.
1 /*
2 ************************************************************
3 * Copyright 2006,2007 by ACcENULLNULL MNRF *
4 * *
5 * http://www.access.edu.au *
6 * Primary Business: Queensland, Australia *
7 * Licensed under the Open NULLoftware License version 3.0 *
8 * http://www.opensource.org/licenses/osl-3.0.php *
9 * *
10 ************************************************************
11 */
12
13
14 /**************************************************************/
15
16 /* assembles the mass matrix in lumped form */
17
18 /* The coefficient D has to be defined on the integration points or not present. */
19
20 /* lumpedMat has to be initialized before the routine is called. */
21
22 /**************************************************************/
23
24 /* Author: gross@access.edu.au */
25 /* Version: $Id:$ */
26
27 /**************************************************************/
28
29 #include "Assemble.h"
30 #include "Util.h"
31 #ifdef _OPENMP
32 #include <omp.h>
33 #endif
34
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 double time0;
45 dim_t dimensions[ESCRIPT_MAX_DATA_RANK], k, e, len_EM_lumpedMat, q, s;
46 type_t funcspace;
47 index_t color,*row_index=NULL;
48 double *S=NULL, *EM_lumpedMat=NULL, *Vol=NULL, *D_p=NULL, *lumpedMat_p=NULL;
49 register double rtmp;
50 size_t len_EM_lumpedMat_size;
51
52 Finley_resetError();
53
54 if (nodes==NULL || elements==NULL) return;
55 if (isEmpty(lumpedMat) || isEmpty(D)) return;
56 if (isEmpty(lumpedMat) && !isEmpty(D)) {
57 Finley_setError(TYPE_ERROR,"Assemble_LumpedSystem: coefficients are non-zero but no lumped matrix is given.");
58 return;
59 }
60 funcspace=getFunctionSpaceType(D);
61 /* check if all function spaces are the same */
62 if (funcspace==FINLEY_ELEMENTS) {
63 reducedIntegrationOrder=FALSE;
64 } else if (funcspace==FINLEY_FACE_ELEMENTS) {
65 reducedIntegrationOrder=FALSE;
66 } else if (funcspace==FINLEY_REDUCED_ELEMENTS) {
67 reducedIntegrationOrder=TRUE;
68 } else if (funcspace==FINLEY_REDUCED_FACE_ELEMENTS) {
69 reducedIntegrationOrder=TRUE;
70 } else {
71 Finley_setError(TYPE_ERROR,"Assemble_LumpedSystem: assemblage failed because of illegal function space.");
72 }
73 if (! Finley_noError()) return;
74
75 /* set all parameters in p*/
76 Assemble_getAssembleParameters(nodes,elements,NULL,lumpedMat, reducedIntegrationOrder, &p);
77 if (! Finley_noError()) return;
78
79 /* check if all function spaces are the same */
80
81 if (! numSamplesEqual(D,p.numQuad,elements->numElements) ) {
82 sprintf(error_msg,"Assemble_LumpedSystem: sample points of coefficient D don't match (%d,%d)",p.numQuad,elements->numElements);
83 Finley_setError(TYPE_ERROR,error_msg);
84 }
85
86 /* check the dimensions: */
87
88 if (p.numEqu==1) {
89 if (!isEmpty(D)) {
90 if (!isDataPointShapeEqual(D,0,dimensions)) {
91 Finley_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 Finley_setError(TYPE_ERROR,error_msg);
101 }
102 }
103 }
104
105 if (Finley_noError()) {
106 lumpedMat_p=getSampleData(lumpedMat,0);
107 len_EM_lumpedMat=p.row_NN*p.numEqu;
108 len_EM_lumpedMat_size=len_EM_lumpedMat*sizeof(double);
109 expandedD=isExpanded(D);
110 S=p.row_jac->ReferenceElement->S;
111
112 #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp)
113 {
114 EM_lumpedMat=THREAD_MEMALLOC(len_EM_lumpedMat,double);
115 row_index=THREAD_MEMALLOC(p.row_NN,index_t);
116 if ( !Finley_checkPtr(EM_lumpedMat) && !Finley_checkPtr(row_index) ) {
117 if (p.numEqu == 1) {
118 if (expandedD) {
119 #ifndef PASO_MPI
120 for (color=elements->minColor;color<=elements->maxColor;color++) {
121 /* open loop over all elements: */
122 #pragma omp for private(e) schedule(static)
123 for(e=0;e<elements->numElements;e++){
124 if (elements->Color[e]==color) {
125 #else
126 {
127 for(e=0;e<elements->numElements;e++) {
128 {
129 #endif
130 Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
131 memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
132 D_p=getSampleData(D,e);
133 for (s=0;s<p.row_NS;s++) {
134 rtmp=0;
135 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[q];
136 EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
137 }
138 for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
139 Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
140 } /* end color check */
141 } /* end element loop */
142 } /* end color loop */
143 } else {
144 #ifndef PASO_MPI
145 for (color=elements->minColor;color<=elements->maxColor;color++) {
146 /* open loop over all elements: */
147 #pragma omp for private(e) schedule(static)
148 for(e=0;e<elements->numElements;e++){
149 if (elements->Color[e]==color) {
150 #else
151 {
152 for(e=0;e<elements->numElements;e++) {
153 {
154 #endif
155 Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
156 memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
157 D_p=getSampleData(D,e);
158 for (s=0;s<p.row_NS;s++) {
159 rtmp=0;
160 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
161 EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp*D_p[0];
162 }
163 for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
164 Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
165 } /* end color check */
166 } /* end element loop */
167 } /* end color loop */
168 }
169 } else {
170 if (expandedD) {
171 #ifndef PASO_MPI
172 for (color=elements->minColor;color<=elements->maxColor;color++) {
173 /* open loop over all elements: */
174 #pragma omp for private(e) schedule(static)
175 for(e=0;e<elements->numElements;e++){
176 if (elements->Color[e]==color) {
177 #else
178 {
179 for(e=0;e<elements->numElements;e++) {
180 {
181 #endif
182 Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
183 memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
184 D_p=getSampleData(D,e);
185 for (s=0;s<p.row_NS;s++) {
186 for (k=0;k<p.numEqu;k++) {
187 rtmp=0.;
188 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[INDEX2(k,q,p.numEqu)];
189 EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
190 }
191 }
192 for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
193 Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
194 } /* end color check */
195 } /* end element loop */
196 } /* end color loop */
197 } else {
198 #ifndef PASO_MPI
199 for (color=elements->minColor;color<=elements->maxColor;color++) {
200 /* open loop over all elements: */
201 #pragma omp for private(e) schedule(static)
202 for(e=0;e<elements->numElements;e++){
203 if (elements->Color[e]==color) {
204 #else
205 {
206 for(e=0;e<elements->numElements;e++) {
207 {
208 #endif
209 Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
210 memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
211 D_p=getSampleData(D,e);
212 for (s=0;s<p.row_NS;s++) {
213 rtmp=0;
214 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
215 for (k=0;k<p.numEqu;k++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp*D_p[k];
216 }
217 for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
218 Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
219 } /* end color check */
220 } /* end element loop */
221 } /* end color loop */
222 }
223 }
224 } /* end of pointer check */
225 THREAD_MEMFREE(EM_lumpedMat);
226 THREAD_MEMFREE(row_index);
227 } /* end parallel region */
228 }
229 #ifdef Finley_TRACE
230 printf("timing: assemblage lumped PDE: %.4e sec\n",Finley_timer()-time0);
231 #endif
232 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26