/[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 100 - (hide annotations)
Wed Dec 15 03:48:48 2004 UTC (14 years, 11 months ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Assemble_PDE.c
File MIME type: text/plain
File size: 15344 byte(s)
*** empty log message ***

1 jgs 82 /* $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 jgs 100 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 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     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 jgs 100 * Revision 1.3 2004/12/15 03:48:44 jgs
346 jgs 97 * *** empty log message ***
347 jgs 82 *
348 jgs 97 * Revision 1.1.1.1 2004/10/26 06:53:57 jgs
349     * initial import of project esys2
350     *
351 jgs 82 * Revision 1.3 2004/07/30 04:37:06 gross
352     * escript and finley are linking now and RecMeshTest.py has been passed
353     *
354     * Revision 1.2 2004/07/21 05:00:54 gross
355     * name changes in DataC
356     *
357     * Revision 1.1 2004/07/02 04:21:13 gross
358     * Finley C code has been included
359     *
360     *
361     */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26