/[escript]/trunk/finley/src/Assemble_LumpedSystem.c
ViewVC logotype

Annotation of /trunk/finley/src/Assemble_LumpedSystem.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1204 - (hide annotations)
Sat Jun 23 11:43:12 2007 UTC (11 years, 10 months ago) by gross
File MIME type: text/plain
File size: 10214 byte(s)
a frame for an improved lumping procedure added.
1 jgs 150 /*
2 elspeth 616 ************************************************************
3 gross 1204 * Copyright 2006,2007 by ACcENULLNULL MNRF *
4 elspeth 616 * *
5     * http://www.access.edu.au *
6     * Primary Business: Queensland, Australia *
7 gross 1204 * Licensed under the Open NULLoftware License version 3.0 *
8 elspeth 616 * http://www.opensource.org/licenses/osl-3.0.php *
9     * *
10     ************************************************************
11 jgs 150 */
12 jgs 82
13 jgs 150
14 jgs 82 /**************************************************************/
15    
16 gross 1204 /* assembles the mass matrix in lumped form */
17 jgs 82
18 gross 1204 /* The coefficient D has to be defined on the integration points or not present. */
19 jgs 82
20 gross 1204 /* lumpedMat has to be initialized before the routine is called. */
21 jgs 82
22     /**************************************************************/
23    
24 jgs 150 /* Author: gross@access.edu.au */
25 gross 1204 /* Version: $Id:$ */
26 jgs 82
27     /**************************************************************/
28    
29     #include "Assemble.h"
30     #include "Util.h"
31     #ifdef _OPENMP
32     #include <omp.h>
33     #endif
34    
35    
36     /**************************************************************/
37    
38 gross 1204 void Finley_Assemble_LumpedSystem(Finley_NodeFile* nodes,Finley_ElementFile* elements, escriptDataC* lumpedMat, escriptDataC* D)
39     {
40 jgs 82
41 gross 1204 bool_t reducedIntegrationOrder=FALSE, expandedD;
42 jgs 150 char error_msg[LenErrorMsg_MAX];
43 gross 798 Assemble_Parameters p;
44 jgs 82 double time0;
45 gross 1204 dim_t dimensions[ESCRIPT_MAX_DATA_RANK], k, e, len_EM_lumpedMat, q, s;
46 gross 1028 type_t funcspace;
47 gross 1204 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 gross 798
52 jgs 150 Finley_resetError();
53 jgs 82
54     if (nodes==NULL || elements==NULL) return;
55 gross 1204 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 gross 798 }
60 gross 1204 funcspace=getFunctionSpaceType(D);
61 jgs 82 /* check if all function spaces are the same */
62 gross 798 if (funcspace==FINLEY_ELEMENTS) {
63     reducedIntegrationOrder=FALSE;
64     } else if (funcspace==FINLEY_FACE_ELEMENTS) {
65     reducedIntegrationOrder=FALSE;
66 gross 1072 } else if (funcspace==FINLEY_REDUCED_ELEMENTS) {
67     reducedIntegrationOrder=TRUE;
68     } else if (funcspace==FINLEY_REDUCED_FACE_ELEMENTS) {
69     reducedIntegrationOrder=TRUE;
70 gross 798 } else {
71 gross 1204 Finley_setError(TYPE_ERROR,"Assemble_LumpedSystem: assemblage failed because of illegal function space.");
72 gross 798 }
73     if (! Finley_noError()) return;
74 jgs 82
75 gross 798 /* set all parameters in p*/
76 gross 1204 Assemble_getAssembleParameters(nodes,elements,NULL,lumpedMat, reducedIntegrationOrder, &p);
77 gross 798 if (! Finley_noError()) return;
78    
79     /* check if all function spaces are the same */
80    
81 jgs 82 if (! numSamplesEqual(D,p.numQuad,elements->numElements) ) {
82 gross 1204 sprintf(error_msg,"Assemble_LumpedSystem: sample points of coefficient D don't match (%d,%d)",p.numQuad,elements->numElements);
83 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
84 jgs 82 }
85    
86     /* check the dimensions: */
87    
88 gross 1204 if (p.numEqu==1) {
89 jgs 82 if (!isEmpty(D)) {
90     if (!isDataPointShapeEqual(D,0,dimensions)) {
91 gross 1204 Finley_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 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
101 jgs 82 }
102     }
103     }
104 gross 1204
105 jgs 150 if (Finley_noError()) {
106 gross 1204 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 gross 798 } else {
170 gross 1204 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 gross 798 }
224 gross 1204 } /* end of pointer check */
225     THREAD_MEMFREE(EM_lumpedMat);
226     THREAD_MEMFREE(row_index);
227     } /* end parallel region */
228 jgs 82 }
229 gross 1204 #ifdef Finley_TRACE
230     printf("timing: assemblage lumped PDE: %.4e sec\n",Finley_timer()-time0);
231     #endif
232 jgs 82 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26