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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 471 - (show annotations)
Fri Jan 27 01:33:02 2006 UTC (13 years, 5 months ago) by jgs
File MIME type: text/plain
File size: 10311 byte(s)
reorganise finley src tree to remove inc dir and src/finley directory

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 /* 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 /* Author: gross@access.edu.au */
35 /* Version: $Id$ */
36
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 char error_msg[LenErrorMsg_MAX];
51 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 Finley_resetError();
58
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 if (! Finley_noError()) return;
65
66 /* this function assumes NS=NN */
67 if (p.NN!=p.NS) {
68 Finley_setError(SYSTEM_ERROR,"__FILE__: for Finley_Assemble_PDE numNodes and numShapes have to be identical.");
69 return;
70 }
71 if (p.numDim!=p.numElementDim) {
72 Finley_setError(SYSTEM_ERROR,"__FILE__: Finley_Assemble_PDE accepts volume elements only.");
73 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 Finley_setError(TYPE_ERROR,"__FILE__: unexpected function space type for coefficient X");
85 }
86 if (! functionSpaceTypeEqual(funcspace,Y) ) {
87 Finley_setError(TYPE_ERROR,"__FILE__: unexpected function space type for coefficient Y");
88 }
89
90 /* check if all function spaces are the same */
91
92 if (! numSamplesEqual(X,p.numQuad,elements->numElements) ) {
93 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 }
96
97 if (! numSamplesEqual(Y,p.numQuad,elements->numElements) ) {
98 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 }
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 sprintf(error_msg,"__FILE__: coefficient X, expected shape (%d,",dimensions[0]);
109 Finley_setError(TYPE_ERROR,error_msg);
110 }
111 }
112 if (!isEmpty(Y)) {
113 if (!isDataPointShapeEqual(Y,0,dimensions)) {
114 Finley_setError(TYPE_ERROR,"__FILE__: coefficient Y, rank 0 expected.");
115 }
116 }
117 } else {
118 if (!isEmpty(X)) {
119 dimensions[0]=p.numEqu;
120 dimensions[1]=p.numDim;
121 if (!isDataPointShapeEqual(X,2,dimensions)) {
122 sprintf(error_msg,"__FILE__: coefficient X, expected shape (%d,%d)",dimensions[0],dimensions[1]);
123 Finley_setError(TYPE_ERROR,error_msg);
124 }
125 }
126 if (!isEmpty(Y)) {
127 dimensions[0]=p.numEqu;
128 if (!isDataPointShapeEqual(Y,1,dimensions)) {
129 sprintf(error_msg,"__FILE__: coefficient Y, expected shape (%d,)",dimensions[0]);
130 Finley_setError(TYPE_ERROR,error_msg);
131 }
132 }
133 }
134
135 if (Finley_noError()) {
136 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 #ifdef Finley_TRACE
219 printf("timing: assemblage PDE Right hand Side: %.4e sec\n",Finley_timer()-time0);
220 #endif
221 }
222 }
223 /*
224 * $Log$
225 * 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 * 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 * 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 * 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 * 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 * 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