/[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 82 - (show annotations)
Tue Oct 26 06:53:54 2004 UTC (14 years, 5 months ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Assemble_PDE.c
File MIME type: text/plain
File size: 15251 byte(s)
Initial revision

1 /* $Id$ */
2
3 /**************************************************************/
4
5 /* assembles the system of numEq PDEs into the stiffness matrix S and right hand side F */
6
7 /* -div(A*grad u)-div(B*u)+C*grad u + D*u= -div X + Y */
8
9 /* -(A_{k,i,m,j} u_m,j)_i-(B_{k,i,m} u_m)_i+C_{k,m,j} u_m,j-D_{k,m} u_m = -(X_{k,i})_i + Y_k */
10
11 /* u has numComp components. */
12
13 /* Shape of the coefficients: */
14
15 /* A = numEqu x numDim x numComp x numDim */
16 /* B = numDim x numEqu x numComp */
17 /* C = numEqu x numDim x numComp */
18 /* D = numEqu x numComp */
19 /* X = numEqu x numDim */
20 /* Y = numEqu */
21
22 /* The coefficients A,B,C,D,X and Y have to be defined on the integartion points or not present (=NULL). */
23
24 /* S and F have to be initialized before the routine is called. S or F can be NULL. In this case the left or */
25 /* the right hand side of the PDE is not processed. */
26
27 /* The routine does not consider any boundary conditions. */
28
29 /**************************************************************/
30
31 /* Copyrights by ACcESS Australia, 2003,2004 */
32 /* author: gross@access.edu.au */
33 /* Version: $Id$ */
34
35 /**************************************************************/
36
37 #include "escript/Data/DataC.h"
38 #include "Finley.h"
39 #include "Assemble.h"
40 #include "NodeFile.h"
41 #include "ElementFile.h"
42 #include "Util.h"
43 #ifdef _OPENMP
44 #include <omp.h>
45 #endif
46
47
48 /**************************************************************/
49
50 void Finley_Assemble_PDE(Finley_NodeFile* nodes,Finley_ElementFile* elements,Finley_SystemMatrix* S, escriptDataC* F,
51 escriptDataC* A, escriptDataC* B, escriptDataC* C, escriptDataC* D, escriptDataC* X, escriptDataC* Y ) {
52
53 double *EM_S=NULL,*EM_F=NULL,*V=NULL,*dVdv=NULL,*dSdV=NULL,*Vol=NULL,*dvdV=NULL;
54 double time0;
55 int dimensions[ESCRIPT_MAX_DATA_RANK],e,q,color;
56 Assemble_Parameters p;
57 maybelong *index_row=NULL,*index_col=NULL;
58
59 if (nodes==NULL || elements==NULL) return;
60 if (S==NULL && isEmpty(F)) return;
61
62 /* set all parameters in p*/
63 Assemble_getAssembleParameters(nodes,elements,S,F,&p);
64 if (Finley_ErrorCode!=NO_ERROR) return;
65
66 /* this function assumes NS=NN */
67 if (p.NN!=p.NS) {
68 Finley_ErrorCode=SYSTEM_ERROR;
69 sprintf(Finley_ErrorMsg,"for Finley_Assemble_PDE numNodes and numShapes have to be identical.");
70 return;
71 }
72 if (p.numDim!=p.numElementDim) {
73 Finley_ErrorCode=SYSTEM_ERROR;
74 sprintf(Finley_ErrorMsg,"Finley_Assemble_PDE accepts volume elements only.");
75 return;
76 }
77 /* get a functionspace */
78 int funcspace=UNKNOWN;
79 updateFunctionSpaceType(funcspace,A);
80 updateFunctionSpaceType(funcspace,B);
81 updateFunctionSpaceType(funcspace,C);
82 updateFunctionSpaceType(funcspace,D);
83 updateFunctionSpaceType(funcspace,X);
84 updateFunctionSpaceType(funcspace,Y);
85 if (funcspace==UNKNOWN) return; /* all data are empty */
86
87 /* check if all function spaces are the same */
88
89 if (! functionSpaceTypeEqual(funcspace,A) ) {
90 Finley_ErrorCode=TYPE_ERROR;
91 sprintf(Finley_ErrorMsg,"unexpected function space typ for coefficient A");
92 }
93 if (! functionSpaceTypeEqual(funcspace,B) ) {
94 Finley_ErrorCode=TYPE_ERROR;
95 sprintf(Finley_ErrorMsg,"unexpected function space typ for coefficient B");
96 }
97 if (! functionSpaceTypeEqual(funcspace,C) ) {
98 Finley_ErrorCode=TYPE_ERROR;
99 sprintf(Finley_ErrorMsg,"unexpected function space typ for coefficient C");
100 }
101 if (! functionSpaceTypeEqual(funcspace,D) ) {
102 Finley_ErrorCode=TYPE_ERROR;
103 sprintf(Finley_ErrorMsg,"unexpected function space typ for coefficient D");
104 }
105 if (! functionSpaceTypeEqual(funcspace,X) ) {
106 Finley_ErrorCode=TYPE_ERROR;
107 sprintf(Finley_ErrorMsg,"unexpected function space typ for coefficient X");
108 }
109 if (! functionSpaceTypeEqual(funcspace,Y) ) {
110 Finley_ErrorCode=TYPE_ERROR;
111 sprintf(Finley_ErrorMsg,"unexpected function space typ for coefficient Y");
112 }
113
114 /* check if all function spaces are the same */
115
116 if (! numSamplesEqual(A,p.numQuad,elements->numElements) ) {
117 Finley_ErrorCode=TYPE_ERROR;
118 sprintf(Finley_ErrorMsg,"sample points of coefficient A don't match (%d,%d)",p.numQuad,elements->numElements);
119 }
120
121 if (! numSamplesEqual(B,p.numQuad,elements->numElements) ) {
122 Finley_ErrorCode=TYPE_ERROR;
123 sprintf(Finley_ErrorMsg,"sample points of coefficient B don't match (%d,%d)",p.numQuad,elements->numElements);
124 }
125
126 if (! numSamplesEqual(C,p.numQuad,elements->numElements) ) {
127 Finley_ErrorCode=TYPE_ERROR;
128 sprintf(Finley_ErrorMsg,"sample points of coefficient C don't match (%d,%d)",p.numQuad,elements->numElements);
129 }
130
131 if (! numSamplesEqual(D,p.numQuad,elements->numElements) ) {
132 Finley_ErrorCode=TYPE_ERROR;
133 sprintf(Finley_ErrorMsg,"sample points of coefficient D don't match (%d,%d)",p.numQuad,elements->numElements);
134 }
135
136 if (! numSamplesEqual(X,p.numQuad,elements->numElements) ) {
137 Finley_ErrorCode=TYPE_ERROR;
138 sprintf(Finley_ErrorMsg,"sample points of coefficient X don't match (%d,%d)",p.numQuad,elements->numElements);
139 }
140
141 if (! numSamplesEqual(Y,p.numQuad,elements->numElements) ) {
142 Finley_ErrorCode=TYPE_ERROR;
143 sprintf(Finley_ErrorMsg,"sample points of coefficient Y don't match (%d,%d)",p.numQuad,elements->numElements);
144 }
145
146 /* check the dimensions: */
147
148 if (p.numEqu==1 && p.numComp==1) {
149 if (!isEmpty(A)) {
150 dimensions[0]=p.numDim;
151 dimensions[1]=p.numDim;
152 if (!isDataPointShapeEqual(A,2,dimensions)) {
153 Finley_ErrorCode=TYPE_ERROR;
154 sprintf(Finley_ErrorMsg,"coefficient A: illegal shape, expected shape (%d,%d)",dimensions[0],dimensions[1]);
155 }
156 }
157 if (!isEmpty(B)) {
158 dimensions[0]=p.numDim;
159 if (!isDataPointShapeEqual(B,1,dimensions)) {
160 Finley_ErrorCode=TYPE_ERROR;
161 sprintf(Finley_ErrorMsg,"coefficient B: illegal shape (%d,)",dimensions[0]);
162 }
163 }
164 if (!isEmpty(C)) {
165 dimensions[0]=p.numDim;
166 if (!isDataPointShapeEqual(C,1,dimensions)) {
167 Finley_ErrorCode=TYPE_ERROR;
168 sprintf(Finley_ErrorMsg,"coefficient C, expected shape (%d,)",dimensions[0]);
169 }
170 }
171 if (!isEmpty(D)) {
172 if (!isDataPointShapeEqual(D,0,dimensions)) {
173 Finley_ErrorCode=TYPE_ERROR;
174 sprintf(Finley_ErrorMsg,"coefficient D, rank 0 expected.");
175 }
176 }
177 if (!isEmpty(X)) {
178 dimensions[0]=p.numDim;
179 if (!isDataPointShapeEqual(X,1,dimensions)) {
180 Finley_ErrorCode=TYPE_ERROR;
181 sprintf(Finley_ErrorMsg,"coefficient X, expected shape (%d,",dimensions[0]);
182 }
183 }
184 if (!isEmpty(Y)) {
185 if (!isDataPointShapeEqual(Y,0,dimensions)) {
186 Finley_ErrorCode=TYPE_ERROR;
187 sprintf(Finley_ErrorMsg,"coefficient Y, rank 0 expected.");
188 }
189 }
190 } else {
191 if (!isEmpty(A)) {
192 dimensions[0]=p.numEqu;
193 dimensions[1]=p.numDim;
194 dimensions[2]=p.numComp;
195 dimensions[3]=p.numDim;
196 if (!isDataPointShapeEqual(A,4,dimensions)) {
197 Finley_ErrorCode=TYPE_ERROR;
198 sprintf(Finley_ErrorMsg,"coefficient A, expected shape (%d,%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2],dimensions[3]);
199 }
200 }
201 if (!isEmpty(B)) {
202 dimensions[0]=p.numEqu;
203 dimensions[1]=p.numDim;
204 dimensions[2]=p.numComp;
205 if (!isDataPointShapeEqual(B,3,dimensions)) {
206 Finley_ErrorCode=TYPE_ERROR;
207 sprintf(Finley_ErrorMsg,"coefficient B, expected shape (%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2]);
208 }
209 }
210 if (!isEmpty(C)) {
211 dimensions[0]=p.numEqu;
212 dimensions[1]=p.numComp;
213 dimensions[2]=p.numDim;
214 if (!isDataPointShapeEqual(C,3,dimensions)) {
215 Finley_ErrorCode=TYPE_ERROR;
216 sprintf(Finley_ErrorMsg,"coefficient C, expected shape (%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2]);
217 }
218 }
219 if (!isEmpty(D)) {
220 dimensions[0]=p.numEqu;
221 dimensions[1]=p.numComp;
222 if (!isDataPointShapeEqual(D,2,dimensions)) {
223 Finley_ErrorCode=TYPE_ERROR;
224 sprintf(Finley_ErrorMsg,"coefficient D, expected shape (%d,%d)",dimensions[0],dimensions[1]);
225 }
226 }
227 if (!isEmpty(X)) {
228 dimensions[0]=p.numEqu;
229 dimensions[1]=p.numDim;
230 if (!isDataPointShapeEqual(X,2,dimensions)) {
231 Finley_ErrorCode=TYPE_ERROR;
232 sprintf(Finley_ErrorMsg,"coefficient X, expected shape (%d,%d)",dimensions[0],dimensions[1]);
233 }
234 }
235 if (!isEmpty(Y)) {
236 dimensions[0]=p.numEqu;
237 if (!isDataPointShapeEqual(Y,1,dimensions)) {
238 Finley_ErrorCode=TYPE_ERROR;
239 sprintf(Finley_ErrorMsg,"coefficient Y, expected shape (%d,)",dimensions[0]);
240 }
241 }
242 }
243
244 if (Finley_ErrorCode==NO_ERROR) {
245 time0=Finley_timer();
246 #pragma omp parallel private(index_col,index_row,EM_S,EM_F,V,dVdv,dSdV,Vol,dvdV,color,q) \
247 firstprivate(elements,nodes,S,F,A,B,C,D,X,Y)
248 {
249 EM_S=EM_F=V=dVdv=dSdV=Vol=dvdV=NULL;
250 index_row=index_col=NULL;
251
252 /* allocate work arrays: */
253
254 EM_S=(double*) THREAD_MEMALLOC(p.NN_row*p.NN_col*p.numEqu*p.numComp*sizeof(double));
255 EM_F=(double*) THREAD_MEMALLOC(p.NN_row*p.numEqu*sizeof(double));
256 V=(double*) THREAD_MEMALLOC(p.NN*p.numDim*sizeof(double));
257 dVdv=(double*) THREAD_MEMALLOC(p.numDim*p.numDim*p.numQuad*sizeof(double));
258 dvdV=(double*) THREAD_MEMALLOC(p.numDim*p.numDim*p.numQuad*sizeof(double));
259 dSdV=(double*) THREAD_MEMALLOC(p.NS_row*p.numQuad*p.numDim*sizeof(double));
260 Vol=(double*) THREAD_MEMALLOC(p.numQuad*sizeof(double));
261 index_col=(maybelong*) THREAD_MEMALLOC(p.NN_col*sizeof(maybelong));
262 index_row=(maybelong*) THREAD_MEMALLOC(p.NN_row*sizeof(maybelong));
263
264 if (! (Finley_checkPtr(EM_S) || Finley_checkPtr(EM_F) || Finley_checkPtr(V) || Finley_checkPtr(index_col) ||
265 Finley_checkPtr(index_row) || Finley_checkPtr(dVdv) || Finley_checkPtr(dSdV) || Finley_checkPtr(Vol) )) {
266
267 /* open loop over all colors: */
268 for (color=0;color<elements->numColors;color++) {
269 /* open loop over all elements: */
270 #pragma omp for private(e) schedule(static)
271 for(e=0;e<elements->numElements;e++){
272 if (elements->Color[e]==color) {
273 for (q=0;q<p.NN_row;q++) index_row[q]=p.label_row[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
274 /* gather V-coordinates of nodes into V: */
275 Finley_Util_Gather_double(p.NN,&(elements->Nodes[INDEX2(0,e,p.NN)]),p.numDim,nodes->Coordinates,V);
276 /* calculate dVdv(i,j,q)=V(i,k)*DSDv(k,j,q) */
277 Finley_Util_SmallMatMult(p.numDim,p.numDim*p.numQuad,dVdv,p.NS,V,p.referenceElement->dSdv);
278 /* dvdV=invert(dVdv) inplace */
279 Finley_Util_InvertSmallMat(p.numQuad,p.numDim,dVdv,dvdV,Vol);
280 /* calculate dSdV=DSDv*DvDV */
281 Finley_Util_SmallMatSetMult(p.numQuad,p.NS_row,p.numDim,dSdV,p.numDim,p.referenceElement_row->dSdv,dvdV);
282 /* scale volume: */
283 for (q=0;q<p.numQuad;q++) Vol[q]=ABS(Vol[q]*p.referenceElement->QuadWeights[q]);
284
285 /* integration for the stiffness matrix: */
286 /* in order to optimze the number of operations the case of constants coefficience needs a bit more work */
287 /* to make use of some symmetry. */
288
289 if (S!=NULL) {
290 for (q=0;q<p.NN_row*p.NN_col*p.numEqu*p.numComp;q++) EM_S[q]=0;
291 if (p.numEqu==1 && p.numComp==1) {
292 Finley_Assemble_PDEMatrix_Single2(p.NS_row,p.numDim,p.numQuad,
293 p.referenceElement_row->S,dSdV,Vol,p.NN_row,EM_S,
294 getSampleData(A,e),isExpanded(A),
295 getSampleData(B,e),isExpanded(B),
296 getSampleData(C,e),isExpanded(C),
297 getSampleData(D,e),isExpanded(D));
298 } else {
299 Finley_Assemble_PDEMatrix_System2(p.NS_row,p.numDim,p.numQuad,p.numEqu,p.numComp,
300 p.referenceElement_row->S,dSdV,Vol,p.NN_row,EM_S,
301 getSampleData(A,e),isExpanded(A),
302 getSampleData(B,e),isExpanded(B),
303 getSampleData(C,e),isExpanded(C),
304 getSampleData(D,e),isExpanded(D));
305 }
306 for (q=0;q<p.NN_col;q++) index_col[q]=p.label_col[elements->Nodes[INDEX2(p.col_node[q],e,p.NN)]];
307 Finley_SystemMatrix_add(S,p.NN_row,index_row,p.numEqu,p.NN_col,index_col,p.numComp,EM_S);
308 }
309 if (!isEmpty(F)) {
310 for (q=0;q<p.NN_row*p.numEqu;q++) EM_F[q]=0;
311 if (p.numEqu==1) {
312 Finley_Assemble_RHSMatrix_Single(p.NS_row,p.numDim,p.numQuad,
313 p.referenceElement_row->S,dSdV,Vol,p.NN_row,EM_F,
314 getSampleData(X,e),isExpanded(X),
315 getSampleData(Y,e),isExpanded(Y));
316 } else {
317 Finley_Assemble_RHSMatrix_System(p.NS_row,p.numDim,p.numQuad,
318 p.numEqu,p.referenceElement_row->S,dSdV,Vol,p.NN_row,EM_F,
319 getSampleData(X,e),isExpanded(X),
320 getSampleData(Y,e),isExpanded(Y));
321 }
322 /* add */
323 Finley_Util_AddScatter(p.NN_row,index_row,p.numEqu,EM_F,getSampleData(F,0));
324 }
325 }
326 }
327 }
328 }
329 /* clean up */
330 THREAD_MEMFREE(EM_S);
331 THREAD_MEMFREE(EM_F);
332 THREAD_MEMFREE(V);
333 THREAD_MEMFREE(dVdv);
334 THREAD_MEMFREE(dvdV);
335 THREAD_MEMFREE(dSdV);
336 THREAD_MEMFREE(Vol);
337 THREAD_MEMFREE(index_col);
338 THREAD_MEMFREE(index_row);
339 }
340 printf("timing: assemblage PDE: %.4e sec\n",Finley_timer()-time0);
341 }
342 }
343 /*
344 * $Log$
345 * Revision 1.1 2004/10/26 06:53:57 jgs
346 * Initial revision
347 *
348 * Revision 1.3 2004/07/30 04:37:06 gross
349 * escript and finley are linking now and RecMeshTest.py has been passed
350 *
351 * Revision 1.2 2004/07/21 05:00:54 gross
352 * name changes in DataC
353 *
354 * Revision 1.1 2004/07/02 04:21:13 gross
355 * Finley C code has been included
356 *
357 *
358 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26