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

Contents of /trunk-mpi-branch/finley/src/Assemble_LumpedSystem.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1253 - (show annotations)
Fri Aug 17 04:09:29 2007 UTC (12 years, 9 months ago) by gross
File MIME type: text/plain
File size: 9535 byte(s)
finally everthing for MPI is in the code. now testing is needed
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: Assemble_LumpedSystem.c 1204 2007-06-23 11:43:12Z gross $ */
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 for (color=elements->minColor;color<=elements->maxColor;color++) {
120 /* open loop over all elements: */
121 #pragma omp for private(e) schedule(static)
122 for(e=0;e<elements->numElements;e++){
123 if (elements->Color[e]==color) {
124 Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
125 memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
126 D_p=getSampleData(D,e);
127 for (s=0;s<p.row_NS;s++) {
128 rtmp=0;
129 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[q];
130 EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
131 }
132 for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
133 Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
134 } /* end color check */
135 } /* end element loop */
136 } /* end color loop */
137 } else {
138 for (color=elements->minColor;color<=elements->maxColor;color++) {
139 /* open loop over all elements: */
140 #pragma omp for private(e) schedule(static)
141 for(e=0;e<elements->numElements;e++){
142 if (elements->Color[e]==color) {
143 Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
144 memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
145 D_p=getSampleData(D,e);
146 for (s=0;s<p.row_NS;s++) {
147 rtmp=0;
148 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
149 EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp*D_p[0];
150 }
151 for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
152 Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
153 } /* end color check */
154 } /* end element loop */
155 } /* end color loop */
156 }
157 } else {
158 if (expandedD) {
159 for (color=elements->minColor;color<=elements->maxColor;color++) {
160 /* open loop over all elements: */
161 #pragma omp for private(e) schedule(static)
162 for(e=0;e<elements->numElements;e++){
163 if (elements->Color[e]==color) {
164 Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
165 memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
166 D_p=getSampleData(D,e);
167 for (s=0;s<p.row_NS;s++) {
168 for (k=0;k<p.numEqu;k++) {
169 rtmp=0.;
170 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[INDEX2(k,q,p.numEqu)];
171 EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
172 }
173 }
174 for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
175 Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
176 } /* end color check */
177 } /* end element loop */
178 } /* end color loop */
179 } else {
180 for (color=elements->minColor;color<=elements->maxColor;color++) {
181 /* open loop over all elements: */
182 #pragma omp for private(e) schedule(static)
183 for(e=0;e<elements->numElements;e++){
184 if (elements->Color[e]==color) {
185 Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
186 memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
187 D_p=getSampleData(D,e);
188 for (s=0;s<p.row_NS;s++) {
189 rtmp=0;
190 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
191 for (k=0;k<p.numEqu;k++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp*D_p[k];
192 }
193 for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
194 Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
195 } /* end color check */
196 } /* end element loop */
197 } /* end color loop */
198 }
199 }
200 } /* end of pointer check */
201 THREAD_MEMFREE(EM_lumpedMat);
202 THREAD_MEMFREE(row_index);
203 } /* end parallel region */
204 }
205 #ifdef Finley_TRACE
206 printf("timing: assemblage lumped PDE: %.4e sec\n",Finley_timer()-time0);
207 #endif
208 }

  ViewVC Help
Powered by ViewVC 1.1.26