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

Contents of /trunk/esys2/finley/src/finleyC/Assemble_PDE_RHS.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 147 - (show annotations)
Fri Aug 12 01:45:47 2005 UTC (14 years, 11 months ago) by jgs
File MIME type: text/plain
File size: 9061 byte(s)
erge of development branch dev-02 back to main trunk on 2005-08-12

1 /* $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 printf("timing: assemblage PDE Right hand Side: %.4e sec\n",Finley_timer()-time0);
215 }
216 }
217 /*
218 * $Log$
219 * Revision 1.2 2005/08/12 01:45:42 jgs
220 * erge of development branch dev-02 back to main trunk on 2005-08-12
221 *
222 * Revision 1.1.2.1 2005/08/04 22:41:11 gross
223 * some extra routines for finley that might speed-up RHS assembling in some cases (not actived right now)
224 *
225 *
226 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26