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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 532 - (hide annotations)
Wed Feb 15 09:45:53 2006 UTC (13 years, 2 months ago) by gross
File MIME type: text/plain
File size: 17175 byte(s)
first steps towards the reuse of the element jacobians
1 jgs 150 /*
2     ******************************************************************************
3     * *
4     * COPYRIGHT ACcESS 2003,2004,2005 - All Rights Reserved *
5     * *
6     * This software is the property of ACcESS. No part of this code *
7     * may be copied in any form or by any means without the expressed written *
8     * consent of ACcESS. Copying, use or modification of this software *
9     * by any unauthorised person is illegal unless that person has a software *
10     * license agreement with ACcESS. *
11     * *
12     ******************************************************************************
13     */
14 jgs 82
15 jgs 150
16 jgs 82 /**************************************************************/
17    
18     /* assembles the system of numEq PDEs into the stiffness matrix S and right hand side F */
19    
20     /* -div(A*grad u)-div(B*u)+C*grad u + D*u= -div X + Y */
21    
22     /* -(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 */
23    
24     /* u has numComp components. */
25    
26     /* Shape of the coefficients: */
27    
28     /* A = numEqu x numDim x numComp x numDim */
29     /* B = numDim x numEqu x numComp */
30     /* C = numEqu x numDim x numComp */
31     /* D = numEqu x numComp */
32     /* X = numEqu x numDim */
33     /* Y = numEqu */
34    
35     /* The coefficients A,B,C,D,X and Y have to be defined on the integartion points or not present (=NULL). */
36    
37     /* S and F have to be initialized before the routine is called. S or F can be NULL. In this case the left or */
38     /* the right hand side of the PDE is not processed. */
39    
40     /* The routine does not consider any boundary conditions. */
41    
42     /**************************************************************/
43    
44 jgs 150 /* Author: gross@access.edu.au */
45     /* Version: $Id$ */
46 jgs 82
47     /**************************************************************/
48    
49     #include "Assemble.h"
50     #include "Util.h"
51     #ifdef _OPENMP
52     #include <omp.h>
53     #endif
54    
55    
56     /**************************************************************/
57    
58 jgs 150 void Finley_Assemble_PDE(Finley_NodeFile* nodes,Finley_ElementFile* elements,Paso_SystemMatrix* S, escriptDataC* F,
59 jgs 82 escriptDataC* A, escriptDataC* B, escriptDataC* C, escriptDataC* D, escriptDataC* X, escriptDataC* Y ) {
60    
61 jgs 150 char error_msg[LenErrorMsg_MAX];
62 jgs 82 double *EM_S=NULL,*EM_F=NULL,*V=NULL,*dVdv=NULL,*dSdV=NULL,*Vol=NULL,*dvdV=NULL;
63     double time0;
64 jgs 123 dim_t dimensions[ESCRIPT_MAX_DATA_RANK],e,q;
65 jgs 82 Assemble_Parameters p;
66 jgs 123 index_t *index_row=NULL,*index_col=NULL,color;
67 jgs 150 Finley_resetError();
68 jgs 82
69     if (nodes==NULL || elements==NULL) return;
70     if (S==NULL && isEmpty(F)) return;
71    
72     /* set all parameters in p*/
73     Assemble_getAssembleParameters(nodes,elements,S,F,&p);
74 jgs 150 if (! Finley_noError()) return;
75 jgs 82
76     /* this function assumes NS=NN */
77     if (p.NN!=p.NS) {
78 gross 532 Finley_setError(SYSTEM_ERROR,"Finley_Assemble_PDE: for Finley_Assemble_PDE numNodes and numShapes have to be identical.");
79 jgs 82 return;
80     }
81     if (p.numDim!=p.numElementDim) {
82 gross 532 Finley_setError(SYSTEM_ERROR,"Finley_Assemble_PDE: Finley_Assemble_PDE accepts volume elements only.");
83 jgs 82 return;
84     }
85     /* get a functionspace */
86 jgs 123 type_t funcspace=UNKNOWN;
87 jgs 82 updateFunctionSpaceType(funcspace,A);
88     updateFunctionSpaceType(funcspace,B);
89     updateFunctionSpaceType(funcspace,C);
90     updateFunctionSpaceType(funcspace,D);
91     updateFunctionSpaceType(funcspace,X);
92     updateFunctionSpaceType(funcspace,Y);
93     if (funcspace==UNKNOWN) return; /* all data are empty */
94    
95     /* check if all function spaces are the same */
96    
97     if (! functionSpaceTypeEqual(funcspace,A) ) {
98 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient A");
99 jgs 82 }
100     if (! functionSpaceTypeEqual(funcspace,B) ) {
101 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient B");
102 jgs 82 }
103     if (! functionSpaceTypeEqual(funcspace,C) ) {
104 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient C");
105 jgs 82 }
106     if (! functionSpaceTypeEqual(funcspace,D) ) {
107 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient D");
108 jgs 82 }
109     if (! functionSpaceTypeEqual(funcspace,X) ) {
110 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient X");
111 jgs 82 }
112     if (! functionSpaceTypeEqual(funcspace,Y) ) {
113 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient Y");
114 jgs 82 }
115    
116     /* check if all function spaces are the same */
117    
118     if (! numSamplesEqual(A,p.numQuad,elements->numElements) ) {
119 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient A don't match (%d,%d)",p.numQuad,elements->numElements);
120 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
121 jgs 82 }
122    
123     if (! numSamplesEqual(B,p.numQuad,elements->numElements) ) {
124 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient B don't match (%d,%d)",p.numQuad,elements->numElements);
125 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
126 jgs 82 }
127    
128     if (! numSamplesEqual(C,p.numQuad,elements->numElements) ) {
129 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient C don't match (%d,%d)",p.numQuad,elements->numElements);
130 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
131 jgs 82 }
132    
133     if (! numSamplesEqual(D,p.numQuad,elements->numElements) ) {
134 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient D don't match (%d,%d)",p.numQuad,elements->numElements);
135 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
136 jgs 82 }
137    
138     if (! numSamplesEqual(X,p.numQuad,elements->numElements) ) {
139 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient X don't match (%d,%d)",p.numQuad,elements->numElements);
140 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
141 jgs 82 }
142    
143     if (! numSamplesEqual(Y,p.numQuad,elements->numElements) ) {
144 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient Y don't match (%d,%d)",p.numQuad,elements->numElements);
145 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
146 jgs 82 }
147    
148     /* check the dimensions: */
149    
150     if (p.numEqu==1 && p.numComp==1) {
151     if (!isEmpty(A)) {
152     dimensions[0]=p.numDim;
153     dimensions[1]=p.numDim;
154     if (!isDataPointShapeEqual(A,2,dimensions)) {
155 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient A: illegal shape, expected shape (%d,%d)",dimensions[0],dimensions[1]);
156 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
157 jgs 82 }
158     }
159     if (!isEmpty(B)) {
160     dimensions[0]=p.numDim;
161     if (!isDataPointShapeEqual(B,1,dimensions)) {
162 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient B: illegal shape (%d,)",dimensions[0]);
163 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
164 jgs 82 }
165     }
166     if (!isEmpty(C)) {
167     dimensions[0]=p.numDim;
168     if (!isDataPointShapeEqual(C,1,dimensions)) {
169 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient C, expected shape (%d,)",dimensions[0]);
170 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
171 jgs 82 }
172     }
173     if (!isEmpty(D)) {
174     if (!isDataPointShapeEqual(D,0,dimensions)) {
175 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficient D, rank 0 expected.");
176 jgs 82 }
177     }
178     if (!isEmpty(X)) {
179     dimensions[0]=p.numDim;
180     if (!isDataPointShapeEqual(X,1,dimensions)) {
181 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient X, expected shape (%d,",dimensions[0]);
182 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
183 jgs 82 }
184     }
185     if (!isEmpty(Y)) {
186     if (!isDataPointShapeEqual(Y,0,dimensions)) {
187 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficient Y, rank 0 expected.");
188 jgs 82 }
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 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient A, expected shape (%d,%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2],dimensions[3]);
198 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
199 jgs 82 }
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 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient B, expected shape (%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2]);
207 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
208 jgs 82 }
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 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient C, expected shape (%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2]);
216 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
217 jgs 82 }
218     }
219     if (!isEmpty(D)) {
220     dimensions[0]=p.numEqu;
221     dimensions[1]=p.numComp;
222     if (!isDataPointShapeEqual(D,2,dimensions)) {
223 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient D, expected shape (%d,%d)",dimensions[0],dimensions[1]);
224 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
225 jgs 82 }
226     }
227     if (!isEmpty(X)) {
228     dimensions[0]=p.numEqu;
229     dimensions[1]=p.numDim;
230     if (!isDataPointShapeEqual(X,2,dimensions)) {
231 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient X, expected shape (%d,%d)",dimensions[0],dimensions[1]);
232 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
233 jgs 82 }
234     }
235     if (!isEmpty(Y)) {
236     dimensions[0]=p.numEqu;
237     if (!isDataPointShapeEqual(Y,1,dimensions)) {
238 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient Y, expected shape (%d,)",dimensions[0]);
239 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
240 jgs 82 }
241     }
242     }
243    
244 jgs 150 if (Finley_noError()) {
245 jgs 82 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 jgs 102 EM_S=(double*) THREAD_MEMALLOC(p.NN_row*p.NN_col*p.numEqu*p.numComp,double);
255     EM_F=(double*) THREAD_MEMALLOC(p.NN_row*p.numEqu,double);
256     V=(double*) THREAD_MEMALLOC(p.NN*p.numDim,double);
257     dVdv=(double*) THREAD_MEMALLOC(p.numDim*p.numDim*p.numQuad,double);
258     dvdV=(double*) THREAD_MEMALLOC(p.numDim*p.numDim*p.numQuad,double);
259     dSdV=(double*) THREAD_MEMALLOC(p.NS_row*p.numQuad*p.numDim,double);
260     Vol=(double*) THREAD_MEMALLOC(p.numQuad,double);
261 jgs 123 index_col=(index_t*) THREAD_MEMALLOC(p.NN_col,index_t);
262     index_row=(index_t*) THREAD_MEMALLOC(p.NN_row,index_t);
263 jgs 82
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 jgs 123 for (color=elements->minColor;color<=elements->maxColor;color++) {
269 jgs 82 /* 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 gross 532 //============================
274 jgs 82 for (q=0;q<p.NN_row;q++) index_row[q]=p.label_row[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
275     /* gather V-coordinates of nodes into V: */
276     Finley_Util_Gather_double(p.NN,&(elements->Nodes[INDEX2(0,e,p.NN)]),p.numDim,nodes->Coordinates,V);
277     /* calculate dVdv(i,j,q)=V(i,k)*DSDv(k,j,q) */
278     Finley_Util_SmallMatMult(p.numDim,p.numDim*p.numQuad,dVdv,p.NS,V,p.referenceElement->dSdv);
279     /* dvdV=invert(dVdv) inplace */
280     Finley_Util_InvertSmallMat(p.numQuad,p.numDim,dVdv,dvdV,Vol);
281     /* calculate dSdV=DSDv*DvDV */
282     Finley_Util_SmallMatSetMult(p.numQuad,p.NS_row,p.numDim,dSdV,p.numDim,p.referenceElement_row->dSdv,dvdV);
283     /* scale volume: */
284     for (q=0;q<p.numQuad;q++) Vol[q]=ABS(Vol[q]*p.referenceElement->QuadWeights[q]);
285 gross 532 //============================
286 jgs 82
287     /* integration for the stiffness matrix: */
288     /* in order to optimze the number of operations the case of constants coefficience needs a bit more work */
289     /* to make use of some symmetry. */
290    
291     if (S!=NULL) {
292     for (q=0;q<p.NN_row*p.NN_col*p.numEqu*p.numComp;q++) EM_S[q]=0;
293     if (p.numEqu==1 && p.numComp==1) {
294     Finley_Assemble_PDEMatrix_Single2(p.NS_row,p.numDim,p.numQuad,
295     p.referenceElement_row->S,dSdV,Vol,p.NN_row,EM_S,
296     getSampleData(A,e),isExpanded(A),
297     getSampleData(B,e),isExpanded(B),
298     getSampleData(C,e),isExpanded(C),
299     getSampleData(D,e),isExpanded(D));
300     } else {
301     Finley_Assemble_PDEMatrix_System2(p.NS_row,p.numDim,p.numQuad,p.numEqu,p.numComp,
302     p.referenceElement_row->S,dSdV,Vol,p.NN_row,EM_S,
303     getSampleData(A,e),isExpanded(A),
304     getSampleData(B,e),isExpanded(B),
305     getSampleData(C,e),isExpanded(C),
306     getSampleData(D,e),isExpanded(D));
307     }
308     for (q=0;q<p.NN_col;q++) index_col[q]=p.label_col[elements->Nodes[INDEX2(p.col_node[q],e,p.NN)]];
309 jgs 150 Finley_Assemble_addToSystemMatrix(S,p.NN_row,index_row,p.numEqu,p.NN_col,index_col,p.numComp,EM_S);
310 jgs 82 }
311     if (!isEmpty(F)) {
312     for (q=0;q<p.NN_row*p.numEqu;q++) EM_F[q]=0;
313     if (p.numEqu==1) {
314     Finley_Assemble_RHSMatrix_Single(p.NS_row,p.numDim,p.numQuad,
315     p.referenceElement_row->S,dSdV,Vol,p.NN_row,EM_F,
316     getSampleData(X,e),isExpanded(X),
317     getSampleData(Y,e),isExpanded(Y));
318     } else {
319     Finley_Assemble_RHSMatrix_System(p.NS_row,p.numDim,p.numQuad,
320     p.numEqu,p.referenceElement_row->S,dSdV,Vol,p.NN_row,EM_F,
321     getSampleData(X,e),isExpanded(X),
322     getSampleData(Y,e),isExpanded(Y));
323     }
324     /* add */
325     Finley_Util_AddScatter(p.NN_row,index_row,p.numEqu,EM_F,getSampleData(F,0));
326     }
327     }
328     }
329     }
330     }
331     /* clean up */
332     THREAD_MEMFREE(EM_S);
333     THREAD_MEMFREE(EM_F);
334     THREAD_MEMFREE(V);
335     THREAD_MEMFREE(dVdv);
336     THREAD_MEMFREE(dvdV);
337     THREAD_MEMFREE(dSdV);
338     THREAD_MEMFREE(Vol);
339     THREAD_MEMFREE(index_col);
340     THREAD_MEMFREE(index_row);
341     }
342 jgs 149 #ifdef Finley_TRACE
343 jgs 82 printf("timing: assemblage PDE: %.4e sec\n",Finley_timer()-time0);
344 jgs 149 #endif
345 jgs 82 }
346     }
347     /*
348     * $Log$
349 jgs 150 * Revision 1.8 2005/09/15 03:44:21 jgs
350     * Merge of development branch dev-02 back to main trunk on 2005-09-15
351     *
352 jgs 149 * Revision 1.7 2005/09/01 03:31:35 jgs
353     * Merge of development branch dev-02 back to main trunk on 2005-09-01
354     *
355 jgs 147 * Revision 1.6 2005/08/12 01:45:42 jgs
356     * erge of development branch dev-02 back to main trunk on 2005-08-12
357     *
358 jgs 150 * Revision 1.5.2.3 2005/09/07 06:26:17 gross
359     * the solver from finley are put into the standalone package paso now
360     *
361 jgs 149 * Revision 1.5.2.2 2005/08/24 02:02:18 gross
362     * timing output switched off. solver output can be swiched through getSolution(verbose=True) now.
363     *
364 jgs 147 * Revision 1.5.2.1 2005/08/03 08:54:27 gross
365     * contact element assemblage was called with wrong element table pointer
366     *
367 jgs 123 * Revision 1.5 2005/07/08 04:07:46 jgs
368     * Merge of development branch back to main trunk on 2005-07-08
369     *
370 jgs 102 * Revision 1.4 2004/12/15 07:08:32 jgs
371 jgs 97 * *** empty log message ***
372 jgs 123 * Revision 1.1.1.1.2.2 2005/06/29 02:34:47 gross
373     * some changes towards 64 integers in finley
374 jgs 82 *
375 jgs 123 * Revision 1.1.1.1.2.1 2004/11/24 01:37:12 gross
376     * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now
377 jgs 97 *
378 jgs 82 *
379 jgs 123 *
380 jgs 82 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26