/[escript]/temp/finley/src/Assemble_jacobeans.c
ViewVC logotype

Contents of /temp/finley/src/Assemble_jacobeans.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 201 - (show annotations)
Wed Nov 23 04:10:21 2005 UTC (13 years, 11 months ago) by jgs
Original Path: trunk/finley/src/finley/Assemble_PDE.c
File MIME type: text/plain
File size: 16827 byte(s)
copy finleyC and CPPAdapter to finley and finley/CPPAdapter to
facilitate scons builds

1 /*
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
15
16 /**************************************************************/
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 /* Author: gross@access.edu.au */
45 /* Version: $Id$ */
46
47 /**************************************************************/
48
49 #include "Assemble.h"
50 #include "Util.h"
51 #ifdef _OPENMP
52 #include <omp.h>
53 #endif
54
55
56 /**************************************************************/
57
58 void Finley_Assemble_PDE(Finley_NodeFile* nodes,Finley_ElementFile* elements,Paso_SystemMatrix* S, escriptDataC* F,
59 escriptDataC* A, escriptDataC* B, escriptDataC* C, escriptDataC* D, escriptDataC* X, escriptDataC* Y ) {
60
61 char error_msg[LenErrorMsg_MAX];
62 double *EM_S=NULL,*EM_F=NULL,*V=NULL,*dVdv=NULL,*dSdV=NULL,*Vol=NULL,*dvdV=NULL;
63 double time0;
64 dim_t dimensions[ESCRIPT_MAX_DATA_RANK],e,q;
65 Assemble_Parameters p;
66 index_t *index_row=NULL,*index_col=NULL,color;
67 Finley_resetError();
68
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 if (! Finley_noError()) return;
75
76 /* this function assumes NS=NN */
77 if (p.NN!=p.NS) {
78 Finley_setError(SYSTEM_ERROR,"__FILE__: for Finley_Assemble_PDE numNodes and numShapes have to be identical.");
79 return;
80 }
81 if (p.numDim!=p.numElementDim) {
82 Finley_setError(SYSTEM_ERROR,"__FILE__: Finley_Assemble_PDE accepts volume elements only.");
83 return;
84 }
85 /* get a functionspace */
86 type_t funcspace=UNKNOWN;
87 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 Finley_setError(TYPE_ERROR,"__FILE__: unexpected function space type for coefficient A");
99 }
100 if (! functionSpaceTypeEqual(funcspace,B) ) {
101 Finley_setError(TYPE_ERROR,"__FILE__: unexpected function space type for coefficient B");
102 }
103 if (! functionSpaceTypeEqual(funcspace,C) ) {
104 Finley_setError(TYPE_ERROR,"__FILE__: unexpected function space type for coefficient C");
105 }
106 if (! functionSpaceTypeEqual(funcspace,D) ) {
107 Finley_setError(TYPE_ERROR,"__FILE__: unexpected function space type for coefficient D");
108 }
109 if (! functionSpaceTypeEqual(funcspace,X) ) {
110 Finley_setError(TYPE_ERROR,"__FILE__: unexpected function space type for coefficient X");
111 }
112 if (! functionSpaceTypeEqual(funcspace,Y) ) {
113 Finley_setError(TYPE_ERROR,"__FILE__: unexpected function space type for coefficient Y");
114 }
115
116 /* check if all function spaces are the same */
117
118 if (! numSamplesEqual(A,p.numQuad,elements->numElements) ) {
119 sprintf(error_msg,"__FILE__: sample points of coefficient A don't match (%d,%d)",p.numQuad,elements->numElements);
120 Finley_setError(TYPE_ERROR,error_msg);
121 }
122
123 if (! numSamplesEqual(B,p.numQuad,elements->numElements) ) {
124 sprintf(error_msg,"__FILE__: sample points of coefficient B don't match (%d,%d)",p.numQuad,elements->numElements);
125 Finley_setError(TYPE_ERROR,error_msg);
126 }
127
128 if (! numSamplesEqual(C,p.numQuad,elements->numElements) ) {
129 sprintf(error_msg,"__FILE__: sample points of coefficient C don't match (%d,%d)",p.numQuad,elements->numElements);
130 Finley_setError(TYPE_ERROR,error_msg);
131 }
132
133 if (! numSamplesEqual(D,p.numQuad,elements->numElements) ) {
134 sprintf(error_msg,"__FILE__: sample points of coefficient D don't match (%d,%d)",p.numQuad,elements->numElements);
135 Finley_setError(TYPE_ERROR,error_msg);
136 }
137
138 if (! numSamplesEqual(X,p.numQuad,elements->numElements) ) {
139 sprintf(error_msg,"__FILE__: sample points of coefficient X don't match (%d,%d)",p.numQuad,elements->numElements);
140 Finley_setError(TYPE_ERROR,error_msg);
141 }
142
143 if (! numSamplesEqual(Y,p.numQuad,elements->numElements) ) {
144 sprintf(error_msg,"__FILE__: sample points of coefficient Y don't match (%d,%d)",p.numQuad,elements->numElements);
145 Finley_setError(TYPE_ERROR,error_msg);
146 }
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 sprintf(error_msg,"__FILE__: coefficient A: illegal shape, expected shape (%d,%d)",dimensions[0],dimensions[1]);
156 Finley_setError(TYPE_ERROR,error_msg);
157 }
158 }
159 if (!isEmpty(B)) {
160 dimensions[0]=p.numDim;
161 if (!isDataPointShapeEqual(B,1,dimensions)) {
162 sprintf(error_msg,"__FILE__: coefficient B: illegal shape (%d,)",dimensions[0]);
163 Finley_setError(TYPE_ERROR,error_msg);
164 }
165 }
166 if (!isEmpty(C)) {
167 dimensions[0]=p.numDim;
168 if (!isDataPointShapeEqual(C,1,dimensions)) {
169 sprintf(error_msg,"__FILE__: coefficient C, expected shape (%d,)",dimensions[0]);
170 Finley_setError(TYPE_ERROR,error_msg);
171 }
172 }
173 if (!isEmpty(D)) {
174 if (!isDataPointShapeEqual(D,0,dimensions)) {
175 Finley_setError(TYPE_ERROR,"__FILE__: coefficient D, rank 0 expected.");
176 }
177 }
178 if (!isEmpty(X)) {
179 dimensions[0]=p.numDim;
180 if (!isDataPointShapeEqual(X,1,dimensions)) {
181 sprintf(error_msg,"__FILE__: coefficient X, expected shape (%d,",dimensions[0]);
182 Finley_setError(TYPE_ERROR,error_msg);
183 }
184 }
185 if (!isEmpty(Y)) {
186 if (!isDataPointShapeEqual(Y,0,dimensions)) {
187 Finley_setError(TYPE_ERROR,"__FILE__: 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 sprintf(error_msg,"__FILE__: coefficient A, expected shape (%d,%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2],dimensions[3]);
198 Finley_setError(TYPE_ERROR,error_msg);
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 sprintf(error_msg,"__FILE__: coefficient B, expected shape (%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2]);
207 Finley_setError(TYPE_ERROR,error_msg);
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 sprintf(error_msg,"__FILE__: coefficient C, expected shape (%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2]);
216 Finley_setError(TYPE_ERROR,error_msg);
217 }
218 }
219 if (!isEmpty(D)) {
220 dimensions[0]=p.numEqu;
221 dimensions[1]=p.numComp;
222 if (!isDataPointShapeEqual(D,2,dimensions)) {
223 sprintf(error_msg,"__FILE__: coefficient D, expected shape (%d,%d)",dimensions[0],dimensions[1]);
224 Finley_setError(TYPE_ERROR,error_msg);
225 }
226 }
227 if (!isEmpty(X)) {
228 dimensions[0]=p.numEqu;
229 dimensions[1]=p.numDim;
230 if (!isDataPointShapeEqual(X,2,dimensions)) {
231 sprintf(error_msg,"__FILE__: coefficient X, expected shape (%d,%d)",dimensions[0],dimensions[1]);
232 Finley_setError(TYPE_ERROR,error_msg);
233 }
234 }
235 if (!isEmpty(Y)) {
236 dimensions[0]=p.numEqu;
237 if (!isDataPointShapeEqual(Y,1,dimensions)) {
238 sprintf(error_msg,"__FILE__: coefficient Y, expected shape (%d,)",dimensions[0]);
239 Finley_setError(TYPE_ERROR,error_msg);
240 }
241 }
242 }
243
244 if (Finley_noError()) {
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,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 index_col=(index_t*) THREAD_MEMALLOC(p.NN_col,index_t);
262 index_row=(index_t*) THREAD_MEMALLOC(p.NN_row,index_t);
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=elements->minColor;color<=elements->maxColor;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_Assemble_addToSystemMatrix(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 #ifdef Finley_TRACE
341 printf("timing: assemblage PDE: %.4e sec\n",Finley_timer()-time0);
342 #endif
343 }
344 }
345 /*
346 * $Log$
347 * Revision 1.8 2005/09/15 03:44:21 jgs
348 * Merge of development branch dev-02 back to main trunk on 2005-09-15
349 *
350 * Revision 1.7 2005/09/01 03:31:35 jgs
351 * Merge of development branch dev-02 back to main trunk on 2005-09-01
352 *
353 * Revision 1.6 2005/08/12 01:45:42 jgs
354 * erge of development branch dev-02 back to main trunk on 2005-08-12
355 *
356 * Revision 1.5.2.3 2005/09/07 06:26:17 gross
357 * the solver from finley are put into the standalone package paso now
358 *
359 * Revision 1.5.2.2 2005/08/24 02:02:18 gross
360 * timing output switched off. solver output can be swiched through getSolution(verbose=True) now.
361 *
362 * Revision 1.5.2.1 2005/08/03 08:54:27 gross
363 * contact element assemblage was called with wrong element table pointer
364 *
365 * Revision 1.5 2005/07/08 04:07:46 jgs
366 * Merge of development branch back to main trunk on 2005-07-08
367 *
368 * Revision 1.4 2004/12/15 07:08:32 jgs
369 * *** empty log message ***
370 * Revision 1.1.1.1.2.2 2005/06/29 02:34:47 gross
371 * some changes towards 64 integers in finley
372 *
373 * Revision 1.1.1.1.2.1 2004/11/24 01:37:12 gross
374 * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now
375 *
376 *
377 *
378 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26