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

Diff of /trunk/finley/src/Assemble_PDE.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 531 by jgs, Fri Jan 27 01:33:02 2006 UTC revision 532 by gross, Wed Feb 15 09:45:53 2006 UTC
# Line 75  void Finley_Assemble_PDE(Finley_NodeFile Line 75  void Finley_Assemble_PDE(Finley_NodeFile
75    
76    /*  this function assumes NS=NN */    /*  this function assumes NS=NN */
77    if (p.NN!=p.NS) {    if (p.NN!=p.NS) {
78      Finley_setError(SYSTEM_ERROR,"__FILE__: for Finley_Assemble_PDE numNodes and numShapes have to be identical.");      Finley_setError(SYSTEM_ERROR,"Finley_Assemble_PDE: for Finley_Assemble_PDE numNodes and numShapes have to be identical.");
79      return;      return;
80    }    }
81    if (p.numDim!=p.numElementDim) {    if (p.numDim!=p.numElementDim) {
82      Finley_setError(SYSTEM_ERROR,"__FILE__: Finley_Assemble_PDE accepts volume elements only.");      Finley_setError(SYSTEM_ERROR,"Finley_Assemble_PDE: Finley_Assemble_PDE accepts volume elements only.");
83      return;      return;
84    }    }
85    /*  get a functionspace */    /*  get a functionspace */
# Line 95  void Finley_Assemble_PDE(Finley_NodeFile Line 95  void Finley_Assemble_PDE(Finley_NodeFile
95    /* check if all function spaces are the same */    /* check if all function spaces are the same */
96    
97    if (! functionSpaceTypeEqual(funcspace,A) ) {    if (! functionSpaceTypeEqual(funcspace,A) ) {
98          Finley_setError(TYPE_ERROR,"__FILE__: unexpected function space type for coefficient A");          Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient A");
99    }    }
100    if (! functionSpaceTypeEqual(funcspace,B) ) {    if (! functionSpaceTypeEqual(funcspace,B) ) {
101          Finley_setError(TYPE_ERROR,"__FILE__: unexpected function space type for coefficient B");          Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient B");
102    }    }
103    if (! functionSpaceTypeEqual(funcspace,C) ) {    if (! functionSpaceTypeEqual(funcspace,C) ) {
104          Finley_setError(TYPE_ERROR,"__FILE__: unexpected function space type for coefficient C");          Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient C");
105    }    }
106    if (! functionSpaceTypeEqual(funcspace,D) ) {    if (! functionSpaceTypeEqual(funcspace,D) ) {
107          Finley_setError(TYPE_ERROR,"__FILE__: unexpected function space type for coefficient D");          Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient D");
108    }    }
109    if (! functionSpaceTypeEqual(funcspace,X) ) {    if (! functionSpaceTypeEqual(funcspace,X) ) {
110          Finley_setError(TYPE_ERROR,"__FILE__: unexpected function space type for coefficient X");          Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient X");
111    }    }
112    if (! functionSpaceTypeEqual(funcspace,Y) ) {    if (! functionSpaceTypeEqual(funcspace,Y) ) {
113          Finley_setError(TYPE_ERROR,"__FILE__: unexpected function space type for coefficient Y");          Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient Y");
114    }    }
115    
116    /* check if all function spaces are the same */    /* check if all function spaces are the same */
117    
118    if (! numSamplesEqual(A,p.numQuad,elements->numElements) ) {    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);          sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient A don't match (%d,%d)",p.numQuad,elements->numElements);
120          Finley_setError(TYPE_ERROR,error_msg);          Finley_setError(TYPE_ERROR,error_msg);
121    }    }
122    
123    if (! numSamplesEqual(B,p.numQuad,elements->numElements) ) {    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);          sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient B don't match (%d,%d)",p.numQuad,elements->numElements);
125          Finley_setError(TYPE_ERROR,error_msg);          Finley_setError(TYPE_ERROR,error_msg);
126    }    }
127    
128    if (! numSamplesEqual(C,p.numQuad,elements->numElements) ) {    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);          sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient C don't match (%d,%d)",p.numQuad,elements->numElements);
130          Finley_setError(TYPE_ERROR,error_msg);          Finley_setError(TYPE_ERROR,error_msg);
131    }    }
132    
133    if (! numSamplesEqual(D,p.numQuad,elements->numElements) ) {    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);          sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient D don't match (%d,%d)",p.numQuad,elements->numElements);
135          Finley_setError(TYPE_ERROR,error_msg);          Finley_setError(TYPE_ERROR,error_msg);
136    }    }
137    
138    if (! numSamplesEqual(X,p.numQuad,elements->numElements) ) {    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);          sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient X don't match (%d,%d)",p.numQuad,elements->numElements);
140          Finley_setError(TYPE_ERROR,error_msg);          Finley_setError(TYPE_ERROR,error_msg);
141    }    }
142    
143    if (! numSamplesEqual(Y,p.numQuad,elements->numElements) ) {    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);          sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient Y don't match (%d,%d)",p.numQuad,elements->numElements);
145          Finley_setError(TYPE_ERROR,error_msg);          Finley_setError(TYPE_ERROR,error_msg);
146    }    }
147    
# Line 152  void Finley_Assemble_PDE(Finley_NodeFile Line 152  void Finley_Assemble_PDE(Finley_NodeFile
152        dimensions[0]=p.numDim;        dimensions[0]=p.numDim;
153        dimensions[1]=p.numDim;        dimensions[1]=p.numDim;
154        if (!isDataPointShapeEqual(A,2,dimensions)) {        if (!isDataPointShapeEqual(A,2,dimensions)) {
155            sprintf(error_msg,"__FILE__: coefficient A: illegal shape, expected shape (%d,%d)",dimensions[0],dimensions[1]);            sprintf(error_msg,"Finley_Assemble_PDE: coefficient A: illegal shape, expected shape (%d,%d)",dimensions[0],dimensions[1]);
156            Finley_setError(TYPE_ERROR,error_msg);            Finley_setError(TYPE_ERROR,error_msg);
157        }        }
158      }      }
159      if (!isEmpty(B)) {      if (!isEmpty(B)) {
160         dimensions[0]=p.numDim;         dimensions[0]=p.numDim;
161         if (!isDataPointShapeEqual(B,1,dimensions)) {         if (!isDataPointShapeEqual(B,1,dimensions)) {
162            sprintf(error_msg,"__FILE__: coefficient B: illegal shape (%d,)",dimensions[0]);            sprintf(error_msg,"Finley_Assemble_PDE: coefficient B: illegal shape (%d,)",dimensions[0]);
163            Finley_setError(TYPE_ERROR,error_msg);            Finley_setError(TYPE_ERROR,error_msg);
164         }         }
165      }      }
166      if (!isEmpty(C)) {      if (!isEmpty(C)) {
167         dimensions[0]=p.numDim;         dimensions[0]=p.numDim;
168         if (!isDataPointShapeEqual(C,1,dimensions)) {         if (!isDataPointShapeEqual(C,1,dimensions)) {
169            sprintf(error_msg,"__FILE__: coefficient C, expected shape (%d,)",dimensions[0]);            sprintf(error_msg,"Finley_Assemble_PDE: coefficient C, expected shape (%d,)",dimensions[0]);
170            Finley_setError(TYPE_ERROR,error_msg);            Finley_setError(TYPE_ERROR,error_msg);
171         }         }
172      }      }
173      if (!isEmpty(D)) {      if (!isEmpty(D)) {
174         if (!isDataPointShapeEqual(D,0,dimensions)) {         if (!isDataPointShapeEqual(D,0,dimensions)) {
175            Finley_setError(TYPE_ERROR,"__FILE__: coefficient D, rank 0 expected.");            Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficient D, rank 0 expected.");
176         }         }
177      }      }
178      if (!isEmpty(X)) {      if (!isEmpty(X)) {
179         dimensions[0]=p.numDim;         dimensions[0]=p.numDim;
180         if (!isDataPointShapeEqual(X,1,dimensions)) {         if (!isDataPointShapeEqual(X,1,dimensions)) {
181            sprintf(error_msg,"__FILE__: coefficient X, expected shape (%d,",dimensions[0]);            sprintf(error_msg,"Finley_Assemble_PDE: coefficient X, expected shape (%d,",dimensions[0]);
182            Finley_setError(TYPE_ERROR,error_msg);            Finley_setError(TYPE_ERROR,error_msg);
183         }         }
184      }      }
185      if (!isEmpty(Y)) {      if (!isEmpty(Y)) {
186         if (!isDataPointShapeEqual(Y,0,dimensions)) {         if (!isDataPointShapeEqual(Y,0,dimensions)) {
187            Finley_setError(TYPE_ERROR,"__FILE__: coefficient Y, rank 0 expected.");            Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficient Y, rank 0 expected.");
188         }         }
189      }      }
190    } else {    } else {
# Line 194  void Finley_Assemble_PDE(Finley_NodeFile Line 194  void Finley_Assemble_PDE(Finley_NodeFile
194        dimensions[2]=p.numComp;        dimensions[2]=p.numComp;
195        dimensions[3]=p.numDim;        dimensions[3]=p.numDim;
196        if (!isDataPointShapeEqual(A,4,dimensions)) {        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]);            sprintf(error_msg,"Finley_Assemble_PDE: coefficient A, expected shape (%d,%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2],dimensions[3]);
198            Finley_setError(TYPE_ERROR,error_msg);            Finley_setError(TYPE_ERROR,error_msg);
199        }        }
200      }      }
# Line 203  void Finley_Assemble_PDE(Finley_NodeFile Line 203  void Finley_Assemble_PDE(Finley_NodeFile
203        dimensions[1]=p.numDim;        dimensions[1]=p.numDim;
204        dimensions[2]=p.numComp;        dimensions[2]=p.numComp;
205        if (!isDataPointShapeEqual(B,3,dimensions)) {        if (!isDataPointShapeEqual(B,3,dimensions)) {
206            sprintf(error_msg,"__FILE__: coefficient B, expected shape (%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2]);            sprintf(error_msg,"Finley_Assemble_PDE: coefficient B, expected shape (%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2]);
207            Finley_setError(TYPE_ERROR,error_msg);            Finley_setError(TYPE_ERROR,error_msg);
208        }        }
209      }      }
# Line 212  void Finley_Assemble_PDE(Finley_NodeFile Line 212  void Finley_Assemble_PDE(Finley_NodeFile
212        dimensions[1]=p.numComp;        dimensions[1]=p.numComp;
213        dimensions[2]=p.numDim;        dimensions[2]=p.numDim;
214        if (!isDataPointShapeEqual(C,3,dimensions)) {        if (!isDataPointShapeEqual(C,3,dimensions)) {
215            sprintf(error_msg,"__FILE__: coefficient C, expected shape (%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2]);            sprintf(error_msg,"Finley_Assemble_PDE: coefficient C, expected shape (%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2]);
216            Finley_setError(TYPE_ERROR,error_msg);            Finley_setError(TYPE_ERROR,error_msg);
217        }        }
218      }      }
# Line 220  void Finley_Assemble_PDE(Finley_NodeFile Line 220  void Finley_Assemble_PDE(Finley_NodeFile
220        dimensions[0]=p.numEqu;        dimensions[0]=p.numEqu;
221        dimensions[1]=p.numComp;        dimensions[1]=p.numComp;
222        if (!isDataPointShapeEqual(D,2,dimensions)) {        if (!isDataPointShapeEqual(D,2,dimensions)) {
223            sprintf(error_msg,"__FILE__: coefficient D, expected shape (%d,%d)",dimensions[0],dimensions[1]);            sprintf(error_msg,"Finley_Assemble_PDE: coefficient D, expected shape (%d,%d)",dimensions[0],dimensions[1]);
224            Finley_setError(TYPE_ERROR,error_msg);            Finley_setError(TYPE_ERROR,error_msg);
225        }        }
226      }      }
# Line 228  void Finley_Assemble_PDE(Finley_NodeFile Line 228  void Finley_Assemble_PDE(Finley_NodeFile
228        dimensions[0]=p.numEqu;        dimensions[0]=p.numEqu;
229        dimensions[1]=p.numDim;        dimensions[1]=p.numDim;
230        if (!isDataPointShapeEqual(X,2,dimensions)) {        if (!isDataPointShapeEqual(X,2,dimensions)) {
231            sprintf(error_msg,"__FILE__: coefficient X, expected shape (%d,%d)",dimensions[0],dimensions[1]);            sprintf(error_msg,"Finley_Assemble_PDE: coefficient X, expected shape (%d,%d)",dimensions[0],dimensions[1]);
232            Finley_setError(TYPE_ERROR,error_msg);            Finley_setError(TYPE_ERROR,error_msg);
233        }        }
234      }      }
235      if (!isEmpty(Y)) {      if (!isEmpty(Y)) {
236        dimensions[0]=p.numEqu;        dimensions[0]=p.numEqu;
237        if (!isDataPointShapeEqual(Y,1,dimensions)) {        if (!isDataPointShapeEqual(Y,1,dimensions)) {
238            sprintf(error_msg,"__FILE__: coefficient Y, expected shape (%d,)",dimensions[0]);            sprintf(error_msg,"Finley_Assemble_PDE: coefficient Y, expected shape (%d,)",dimensions[0]);
239            Finley_setError(TYPE_ERROR,error_msg);            Finley_setError(TYPE_ERROR,error_msg);
240        }        }
241      }      }
# Line 270  void Finley_Assemble_PDE(Finley_NodeFile Line 270  void Finley_Assemble_PDE(Finley_NodeFile
270                #pragma omp for private(e) schedule(static)                #pragma omp for private(e) schedule(static)
271                for(e=0;e<elements->numElements;e++){                for(e=0;e<elements->numElements;e++){
272                  if (elements->Color[e]==color) {                  if (elements->Color[e]==color) {
273    //============================
274                    for (q=0;q<p.NN_row;q++) index_row[q]=p.label_row[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];                    for (q=0;q<p.NN_row;q++) index_row[q]=p.label_row[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
275                    /* gather V-coordinates of nodes into V: */                    /* gather V-coordinates of nodes into V: */
276            Finley_Util_Gather_double(p.NN,&(elements->Nodes[INDEX2(0,e,p.NN)]),p.numDim,nodes->Coordinates,V);            Finley_Util_Gather_double(p.NN,&(elements->Nodes[INDEX2(0,e,p.NN)]),p.numDim,nodes->Coordinates,V);
# Line 281  void Finley_Assemble_PDE(Finley_NodeFile Line 282  void Finley_Assemble_PDE(Finley_NodeFile
282            Finley_Util_SmallMatSetMult(p.numQuad,p.NS_row,p.numDim,dSdV,p.numDim,p.referenceElement_row->dSdv,dvdV);            Finley_Util_SmallMatSetMult(p.numQuad,p.NS_row,p.numDim,dSdV,p.numDim,p.referenceElement_row->dSdv,dvdV);
283                    /*  scale volume: */                    /*  scale volume: */
284            for (q=0;q<p.numQuad;q++) Vol[q]=ABS(Vol[q]*p.referenceElement->QuadWeights[q]);            for (q=0;q<p.numQuad;q++) Vol[q]=ABS(Vol[q]*p.referenceElement->QuadWeights[q]);
285    //============================
286            
287                     /*   integration for the stiffness matrix: */                     /*   integration for the stiffness matrix: */
288                     /*   in order to optimze the number of operations the case of constants coefficience needs a bit more work */                     /*   in order to optimze the number of operations the case of constants coefficience needs a bit more work */

Legend:
Removed from v.531  
changed lines
  Added in v.532

  ViewVC Help
Powered by ViewVC 1.1.26