/[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 149 - (hide annotations)
Thu Sep 1 03:31:39 2005 UTC (14 years, 10 months ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Assemble_PDE_RHS.c
File MIME type: text/plain
File size: 9364 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-01

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26