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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 201 - (hide annotations)
Wed Nov 23 04:10:21 2005 UTC (14 years ago) by jgs
File MIME type: text/plain
File size: 10311 byte(s)
copy finleyC and CPPAdapter to finley and finley/CPPAdapter to
facilitate scons builds

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 147
15     /**************************************************************/
16    
17     /* assembles the system of numEq PDEs into side F */
18    
19     /* -div X + Y */
20    
21     /* Shape of the coefficients: */
22    
23     /* X = numEqu x numDim */
24     /* Y = numEqu */
25    
26     /* The coefficients X and Y have to be defined on the integartion points or not present (=NULL). */
27    
28     /* F has to be initialized before the routine is called. F can be NULL. */
29    
30     /* The routine does not consider any boundary conditions. */
31    
32     /**************************************************************/
33    
34 jgs 150 /* Author: gross@access.edu.au */
35     /* Version: $Id$ */
36 jgs 147
37     /**************************************************************/
38    
39     #include "Assemble.h"
40     #include "Util.h"
41     #ifdef _OPENMP
42     #include <omp.h>
43     #endif
44    
45    
46     /**************************************************************/
47    
48     void Finley_Assemble_PDE_RHS(Finley_NodeFile* nodes,Finley_ElementFile* elements,escriptDataC* F,escriptDataC* X, escriptDataC* Y ) {
49    
50 jgs 150 char error_msg[LenErrorMsg_MAX];
51 jgs 147 double *EM_F=NULL,*V=NULL,*dVdv=NULL,*dSdV=NULL,*Vol=NULL,*dvdV=NULL;
52     double time0;
53     dim_t dimensions[ESCRIPT_MAX_DATA_RANK],e,q;
54     Assemble_Parameters p;
55     index_t *index_row=NULL,color;
56     bool_t assemble;
57 jgs 150 Finley_resetError();
58 jgs 147
59     if (nodes==NULL || elements==NULL) return;
60     if (isEmpty(F) || (isEmpty(Y) && isEmpty(X)) ) return;
61    
62     /* set all parameters in p*/
63     Assemble_getAssembleParameters(nodes,elements,NULL,F,&p);
64 jgs 150 if (! Finley_noError()) return;
65 jgs 147
66     /* this function assumes NS=NN */
67     if (p.NN!=p.NS) {
68 jgs 150 Finley_setError(SYSTEM_ERROR,"__FILE__: for Finley_Assemble_PDE numNodes and numShapes have to be identical.");
69 jgs 147 return;
70     }
71     if (p.numDim!=p.numElementDim) {
72 jgs 150 Finley_setError(SYSTEM_ERROR,"__FILE__: Finley_Assemble_PDE accepts volume elements only.");
73 jgs 147 return;
74     }
75     /* get a functionspace */
76     type_t funcspace=UNKNOWN;
77     updateFunctionSpaceType(funcspace,X);
78     updateFunctionSpaceType(funcspace,Y);
79     if (funcspace==UNKNOWN) return; /* all data are empty */
80    
81     /* check if all function spaces are the same */
82    
83     if (! functionSpaceTypeEqual(funcspace,X) ) {
84 jgs 150 Finley_setError(TYPE_ERROR,"__FILE__: unexpected function space type for coefficient X");
85 jgs 147 }
86     if (! functionSpaceTypeEqual(funcspace,Y) ) {
87 jgs 150 Finley_setError(TYPE_ERROR,"__FILE__: unexpected function space type for coefficient Y");
88 jgs 147 }
89    
90     /* check if all function spaces are the same */
91    
92     if (! numSamplesEqual(X,p.numQuad,elements->numElements) ) {
93 jgs 150 sprintf(error_msg,"__FILE__: sample points of coefficient X don't match (%d,%d)",p.numQuad,elements->numElements);
94     Finley_setError(TYPE_ERROR,error_msg);
95 jgs 147 }
96    
97     if (! numSamplesEqual(Y,p.numQuad,elements->numElements) ) {
98 jgs 150 sprintf(error_msg,"__FILE__: sample points of coefficient Y don't match (%d,%d)",p.numQuad,elements->numElements);
99     Finley_setError(TYPE_ERROR,error_msg);
100 jgs 147 }
101    
102     /* check the dimensions: */
103    
104     if (p.numEqu==1 && p.numComp==1) {
105     if (!isEmpty(X)) {
106     dimensions[0]=p.numDim;
107     if (!isDataPointShapeEqual(X,1,dimensions)) {
108 jgs 150 sprintf(error_msg,"__FILE__: coefficient X, expected shape (%d,",dimensions[0]);
109     Finley_setError(TYPE_ERROR,error_msg);
110 jgs 147 }
111     }
112     if (!isEmpty(Y)) {
113     if (!isDataPointShapeEqual(Y,0,dimensions)) {
114 jgs 150 Finley_setError(TYPE_ERROR,"__FILE__: coefficient Y, rank 0 expected.");
115 jgs 147 }
116     }
117     } else {
118     if (!isEmpty(X)) {
119     dimensions[0]=p.numEqu;
120     dimensions[1]=p.numDim;
121     if (!isDataPointShapeEqual(X,2,dimensions)) {
122 jgs 150 sprintf(error_msg,"__FILE__: coefficient X, expected shape (%d,%d)",dimensions[0],dimensions[1]);
123     Finley_setError(TYPE_ERROR,error_msg);
124 jgs 147 }
125     }
126     if (!isEmpty(Y)) {
127     dimensions[0]=p.numEqu;
128     if (!isDataPointShapeEqual(Y,1,dimensions)) {
129 jgs 150 sprintf(error_msg,"__FILE__: coefficient Y, expected shape (%d,)",dimensions[0]);
130     Finley_setError(TYPE_ERROR,error_msg);
131 jgs 147 }
132     }
133     }
134    
135 jgs 150 if (Finley_noError()) {
136 jgs 147 time0=Finley_timer();
137     #pragma omp parallel private(index_row,EM_F,V,dVdv,dSdV,Vol,dvdV,color,q,assemble)
138     {
139     EM_F=V=dVdv=dSdV=Vol=dvdV=NULL;
140     index_row=NULL;
141    
142     /* allocate work arrays: */
143    
144     EM_F=(double*) THREAD_MEMALLOC(p.NN_row*p.numEqu,double);
145     V=(double*) THREAD_MEMALLOC(p.NN*p.numDim,double);
146     dVdv=(double*) THREAD_MEMALLOC(p.numDim*p.numDim*p.numQuad,double);
147     dvdV=(double*) THREAD_MEMALLOC(p.numDim*p.numDim*p.numQuad,double);
148     dSdV=(double*) THREAD_MEMALLOC(p.NS_row*p.numQuad*p.numDim,double);
149     Vol=(double*) THREAD_MEMALLOC(p.numQuad,double);
150     index_row=(index_t*) THREAD_MEMALLOC(p.NN_row,index_t);
151    
152     if (! (Finley_checkPtr(EM_F) || Finley_checkPtr(V) ||
153     Finley_checkPtr(index_row) || Finley_checkPtr(dVdv) || Finley_checkPtr(dSdV) || Finley_checkPtr(Vol) )) {
154    
155     /* open loop over all colors: */
156     for (color=elements->minColor;color<=elements->maxColor;color++) {
157     /* open loop over all elements: */
158     #pragma omp for private(e) schedule(dynamic)
159     for(e=0;e<elements->numElements;e++){
160     if (elements->Color[e]==color) {
161     /* check if element has to be processed */
162     if (isEmpty(X)) {
163     assemble=FALSE;
164     } else {
165     if (isExpanded(X)) {
166     assemble=Finley_Util_anyNonZeroDouble(p.numEqu*p.numDim*p.numQuad,getSampleData(X,e));
167     } else {
168     assemble=Finley_Util_anyNonZeroDouble(p.numEqu*p.numDim,getSampleData(X,e));
169     }
170     }
171     if (!assemble && !isEmpty(Y)) {
172     if (isExpanded(Y)) {
173     assemble=Finley_Util_anyNonZeroDouble(p.numEqu*p.numQuad,getSampleData(Y,e));
174     } else {
175     assemble=Finley_Util_anyNonZeroDouble(p.numEqu,getSampleData(Y,e));
176     }
177     }
178     if (assemble) {
179     for (q=0;q<p.NN_row;q++) index_row[q]=p.label_row[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
180     /* gather V-coordinates of nodes into V: */
181     Finley_Util_Gather_double(p.NN,&(elements->Nodes[INDEX2(0,e,p.NN)]),p.numDim,nodes->Coordinates,V);
182     /* calculate dVdv(i,j,q)=V(i,k)*DSDv(k,j,q) */
183     Finley_Util_SmallMatMult(p.numDim,p.numDim*p.numQuad,dVdv,p.NS,V,p.referenceElement->dSdv);
184     /* dvdV=invert(dVdv) inplace */
185     Finley_Util_InvertSmallMat(p.numQuad,p.numDim,dVdv,dvdV,Vol);
186     /* calculate dSdV=DSDv*DvDV */
187     Finley_Util_SmallMatSetMult(p.numQuad,p.NS_row,p.numDim,dSdV,p.numDim,p.referenceElement_row->dSdv,dvdV);
188     /* scale volume: */
189     for (q=0;q<p.numQuad;q++) Vol[q]=ABS(Vol[q]*p.referenceElement->QuadWeights[q]);
190     for (q=0;q<p.NN_row*p.numEqu;q++) EM_F[q]=0;
191     if (p.numEqu==1) {
192     Finley_Assemble_RHSMatrix_Single(p.NS_row,p.numDim,p.numQuad,
193     p.referenceElement_row->S,dSdV,Vol,p.NN_row,EM_F,
194     getSampleData(X,e),isExpanded(X),
195     getSampleData(Y,e),isExpanded(Y));
196     } else {
197     Finley_Assemble_RHSMatrix_System(p.NS_row,p.numDim,p.numQuad,
198     p.numEqu,p.referenceElement_row->S,dSdV,Vol,p.NN_row,EM_F,
199     getSampleData(X,e),isExpanded(X),
200     getSampleData(Y,e),isExpanded(Y));
201     }
202     /* add */
203     Finley_Util_AddScatter(p.NN_row,index_row,p.numEqu,EM_F,getSampleData(F,0));
204     }
205     }
206     }
207     }
208     }
209     /* clean up */
210     THREAD_MEMFREE(EM_F);
211     THREAD_MEMFREE(V);
212     THREAD_MEMFREE(dVdv);
213     THREAD_MEMFREE(dvdV);
214     THREAD_MEMFREE(dSdV);
215     THREAD_MEMFREE(Vol);
216     THREAD_MEMFREE(index_row);
217     }
218 jgs 149 #ifdef Finley_TRACE
219 jgs 147 printf("timing: assemblage PDE Right hand Side: %.4e sec\n",Finley_timer()-time0);
220 jgs 149 #endif
221 jgs 147 }
222     }
223     /*
224     * $Log$
225 jgs 150 * Revision 1.4 2005/09/15 03:44:21 jgs
226     * Merge of development branch dev-02 back to main trunk on 2005-09-15
227     *
228 jgs 149 * Revision 1.3 2005/09/01 03:31:35 jgs
229     * Merge of development branch dev-02 back to main trunk on 2005-09-01
230     *
231 jgs 147 * Revision 1.2 2005/08/12 01:45:42 jgs
232     * erge of development branch dev-02 back to main trunk on 2005-08-12
233     *
234 jgs 150 * Revision 1.1.2.3 2005/09/07 06:26:17 gross
235     * the solver from finley are put into the standalone package paso now
236     *
237 jgs 149 * Revision 1.1.2.2 2005/08/24 02:02:18 gross
238     * timing output switched off. solver output can be swiched through getSolution(verbose=True) now.
239     *
240 jgs 147 * Revision 1.1.2.1 2005/08/04 22:41:11 gross
241     * some extra routines for finley that might speed-up RHS assembling in some cases (not actived right now)
242     *
243     *
244     */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26