/[escript]/branches/domexper/dudley/src/Assemble_PDE.c
ViewVC logotype

Diff of /branches/domexper/dudley/src/Assemble_PDE.c

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

revision 971 by ksteube, Wed Feb 14 04:40:49 2007 UTC revision 2748 by gross, Tue Nov 17 07:32:59 2009 UTC
# Line 1  Line 1 
1  /*  
2   ************************************************************  /*******************************************************
3   *          Copyright 2006 by ACcESS MNRF                   *  *
4   *                                                          *  * Copyright (c) 2003-2009 by University of Queensland
5   *              http://www.access.edu.au                    *  * Earth Systems Science Computational Center (ESSCC)
6   *       Primary Business: Queensland, Australia            *  * http://www.uq.edu.au/esscc
7   *  Licensed under the Open Software License version 3.0    *  *
8   *     http://www.opensource.org/licenses/osl-3.0.php       *  * Primary Business: Queensland, Australia
9   *                                                          *  * Licensed under the Open Software License version 3.0
10   ************************************************************  * http://www.opensource.org/licenses/osl-3.0.php
11  */  *
12    *******************************************************/
13    
14    
15  /**************************************************************/  /**************************************************************/
# Line 39  Line 40 
40    
41  /**************************************************************/  /**************************************************************/
42    
 /*  Author: gross@access.edu.au */  
 /*  Version: $Id$ */  
   
 /**************************************************************/  
   
43  #include "Assemble.h"  #include "Assemble.h"
44  #include "Util.h"  #include "Util.h"
45    #include "esysUtils/blocktimer.h"
46  #ifdef _OPENMP  #ifdef _OPENMP
47  #include <omp.h>  #include <omp.h>
48  #endif  #endif
# Line 61  void Finley_Assemble_PDE(Finley_NodeFile Line 58  void Finley_Assemble_PDE(Finley_NodeFile
58    Assemble_Parameters p;    Assemble_Parameters p;
59    double time0;    double time0;
60    dim_t dimensions[ESCRIPT_MAX_DATA_RANK];    dim_t dimensions[ESCRIPT_MAX_DATA_RANK];
61      type_t funcspace;
62      double blocktimer_start = blocktimer_time();
63    
64    Finley_resetError();    Finley_resetError();
65    
66      {
67    #ifdef Paso_MPI
68      int iam, numCPUs;
69      iam = elements->elementDistribution->MPIInfo->rank;
70      numCPUs = elements->elementDistribution->MPIInfo->size;
71    #endif
72      }
73    
74    if (nodes==NULL || elements==NULL) return;    if (nodes==NULL || elements==NULL) return;
75    if (S==NULL && isEmpty(F)) return;    if (S==NULL && isEmpty(F)) return;
76    
77    if (isEmpty(F) && !isEmpty(X) && !isEmpty(F)) {    if (isEmpty(F) && !isEmpty(X) && !isEmpty(F)) {
78          Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: right hand side coefficients are non-zero bat no right hand side vector given.");          Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: right hand side coefficients are non-zero bat no right hand side vector given.");
79    }    }
80    
81    if (S==NULL && !isEmpty(A) && !isEmpty(B) && !isEmpty(C) && !isEmpty(D)) {    if (S==NULL && !isEmpty(A) && !isEmpty(B) && !isEmpty(C) && !isEmpty(D)) {
82          Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficients are non-zero but no matrix is given.");          Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficients are non-zero but no matrix is given.");
83    }    }
84    
85    /*  get the functionspace for this assemblage call */    /*  get the functionspace for this assemblage call */
86    type_t funcspace=UNKNOWN;    funcspace=UNKNOWN;
87    updateFunctionSpaceType(funcspace,A);    updateFunctionSpaceType(funcspace,A);
88    updateFunctionSpaceType(funcspace,B);    updateFunctionSpaceType(funcspace,B);
89    updateFunctionSpaceType(funcspace,C);    updateFunctionSpaceType(funcspace,C);
# Line 85  void Finley_Assemble_PDE(Finley_NodeFile Line 92  void Finley_Assemble_PDE(Finley_NodeFile
92    updateFunctionSpaceType(funcspace,Y);    updateFunctionSpaceType(funcspace,Y);
93    if (funcspace==UNKNOWN) return; /* all  data are empty */    if (funcspace==UNKNOWN) return; /* all  data are empty */
94    
   
95    /* check if all function spaces are the same */    /* check if all function spaces are the same */
96    if (! functionSpaceTypeEqual(funcspace,A) ) {    if (! functionSpaceTypeEqual(funcspace,A) ) {
97          Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient A");          Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient A");
# Line 116  void Finley_Assemble_PDE(Finley_NodeFile Line 122  void Finley_Assemble_PDE(Finley_NodeFile
122         reducedIntegrationOrder=FALSE;         reducedIntegrationOrder=FALSE;
123    } else if (funcspace==FINLEY_CONTACT_ELEMENTS_2)  {    } else if (funcspace==FINLEY_CONTACT_ELEMENTS_2)  {
124         reducedIntegrationOrder=FALSE;         reducedIntegrationOrder=FALSE;
125      } else if (funcspace==FINLEY_REDUCED_ELEMENTS) {
126           reducedIntegrationOrder=TRUE;
127      } else if (funcspace==FINLEY_REDUCED_FACE_ELEMENTS)  {
128           reducedIntegrationOrder=TRUE;
129      } else if (funcspace==FINLEY_REDUCED_CONTACT_ELEMENTS_1)  {
130           reducedIntegrationOrder=TRUE;
131      } else if (funcspace==FINLEY_REDUCED_CONTACT_ELEMENTS_2)  {
132           reducedIntegrationOrder=TRUE;
133    } else {    } else {
134         Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: assemblage failed because of illegal function space.");         Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: assemblage failed because of illegal function space.");
135    }    }
# Line 127  void Finley_Assemble_PDE(Finley_NodeFile Line 141  void Finley_Assemble_PDE(Finley_NodeFile
141    
142    /* check if all function spaces are the same */    /* check if all function spaces are the same */
143    
144    if (! numSamplesEqual(A,p.numQuad,elements->numElements) ) {    if (! numSamplesEqual(A,p.numQuadTotal,elements->numElements) ) {
145          sprintf(error_msg,"Finley_Assemble_PDE: 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.numQuadTotal,elements->numElements);
146          Finley_setError(TYPE_ERROR,error_msg);          Finley_setError(TYPE_ERROR,error_msg);
147    }    }
148    
149    if (! numSamplesEqual(B,p.numQuad,elements->numElements) ) {    if (! numSamplesEqual(B,p.numQuadTotal,elements->numElements) ) {
150          sprintf(error_msg,"Finley_Assemble_PDE: 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.numQuadTotal,elements->numElements);
151          Finley_setError(TYPE_ERROR,error_msg);          Finley_setError(TYPE_ERROR,error_msg);
152    }    }
153    
154    if (! numSamplesEqual(C,p.numQuad,elements->numElements) ) {    if (! numSamplesEqual(C,p.numQuadTotal,elements->numElements) ) {
155          sprintf(error_msg,"Finley_Assemble_PDE: 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.numQuadTotal,elements->numElements);
156          Finley_setError(TYPE_ERROR,error_msg);          Finley_setError(TYPE_ERROR,error_msg);
157    }    }
158    
159    if (! numSamplesEqual(D,p.numQuad,elements->numElements) ) {    if (! numSamplesEqual(D,p.numQuadTotal,elements->numElements) ) {
160          sprintf(error_msg,"Finley_Assemble_PDE: 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.numQuadTotal,elements->numElements);
161          Finley_setError(TYPE_ERROR,error_msg);          Finley_setError(TYPE_ERROR,error_msg);
162    }    }
163    
164    if (! numSamplesEqual(X,p.numQuad,elements->numElements) ) {    if (! numSamplesEqual(X,p.numQuadTotal,elements->numElements) ) {
165          sprintf(error_msg,"Finley_Assemble_PDE: 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.numQuadTotal,elements->numElements);
166          Finley_setError(TYPE_ERROR,error_msg);          Finley_setError(TYPE_ERROR,error_msg);
167    }    }
168    
169    if (! numSamplesEqual(Y,p.numQuad,elements->numElements) ) {    if (! numSamplesEqual(Y,p.numQuadTotal,elements->numElements) ) {
170          sprintf(error_msg,"Finley_Assemble_PDE: 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.numQuadTotal,elements->numElements);
171          Finley_setError(TYPE_ERROR,error_msg);          Finley_setError(TYPE_ERROR,error_msg);
172    }    }
173    
174    /*  check the dimensions: */    /*  check the dimensions: */
175      
176    if (p.numEqu==1 && p.numComp==1) {    if (p.numEqu==1 && p.numComp==1) {
177      if (!isEmpty(A)) {      if (!isEmpty(A)) {
178        dimensions[0]=p.numDim;        dimensions[0]=p.numDim;
# Line 198  void Finley_Assemble_PDE(Finley_NodeFile Line 212  void Finley_Assemble_PDE(Finley_NodeFile
212         if (!isDataPointShapeEqual(Y,0,dimensions)) {         if (!isDataPointShapeEqual(Y,0,dimensions)) {
213            Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficient Y, rank 0 expected.");            Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficient Y, rank 0 expected.");
214         }         }
215      }      }
216    } else {    } else {
217      if (!isEmpty(A)) {      if (!isEmpty(A)) {
218        dimensions[0]=p.numEqu;        dimensions[0]=p.numEqu;
# Line 258  void Finley_Assemble_PDE(Finley_NodeFile Line 272  void Finley_Assemble_PDE(Finley_NodeFile
272          if (p.numEqu > 1) {          if (p.numEqu > 1) {
273            /* system of PDESs */            /* system of PDESs */
274            if (p.numDim==3) {            if (p.numDim==3) {
275              if (p.row_NS == p.col_NS && p.row_NS == p.row_NN && p.col_NS == p.col_NN ) {              if ( p.numSides == 1 ) {
276                 Finley_Assemble_PDE_System2_3D(p,elements,S,F,A,B,C,D,X,Y);                 Finley_Assemble_PDE_System2_3D(p,elements,S,F,A,B,C,D,X,Y);
277              } else if ( p.row_NS == p.col_NS &&  2*p.row_NS == p.row_NN && 2*p.col_NS == p.col_NN ) {              } else if ( p.numSides == 2 ) {
278                 if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {                 if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
279                    Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");                    Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
280                 } else {                 } else {
# Line 270  void Finley_Assemble_PDE(Finley_NodeFile Line 284  void Finley_Assemble_PDE(Finley_NodeFile
284                 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");                 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
285              }              }
286            } else if (p.numDim==2) {            } else if (p.numDim==2) {
287              if ((p.row_NS == p.col_NS) && (p.row_NS == p.row_NN) && (p.col_NS == p.col_NN )) {              if ( p.numSides == 1 ) {
288                 Finley_Assemble_PDE_System2_2D(p,elements,S,F,A,B,C,D,X,Y);                 Finley_Assemble_PDE_System2_2D(p,elements,S,F,A,B,C,D,X,Y);
289              } else if ( p.row_NS == p.col_NS &&  2*p.row_NS == p.row_NN && 2*p.col_NS == p.col_NN ) {              } else if (  p.numSides == 2 ) {
290                 if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {                 if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
291                    Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");                    Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
292                 } else {                 } else {
# Line 281  void Finley_Assemble_PDE(Finley_NodeFile Line 295  void Finley_Assemble_PDE(Finley_NodeFile
295              } else {              } else {
296                 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");                 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
297              }              }
298            } else if (p.numDim==2) {            } else if (p.numDim==1) {
299              if (p.row_NS == p.col_NS && p.row_NS == p.row_NN && p.col_NS == p.col_NN ) {              if ( p.numSides == 1  ) {
300                 Finley_Assemble_PDE_System2_1D(p,elements,S,F,A,B,C,D,X,Y);                 Finley_Assemble_PDE_System2_1D(p,elements,S,F,A,B,C,D,X,Y);
301              } else if ( p.row_NS == p.col_NS &&  2*p.row_NS == p.row_NN && 2*p.col_NS == p.col_NN ) {              } else if ( p.numSides == 2 ) {
302                 if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {                 if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
303                    Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");                    Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
304                 } else {                 } else {
# Line 299  void Finley_Assemble_PDE(Finley_NodeFile Line 313  void Finley_Assemble_PDE(Finley_NodeFile
313          } else {          } else {
314            /* single PDES */            /* single PDES */
315            if (p.numDim==3) {            if (p.numDim==3) {
316              if (p.row_NS == p.col_NS && p.row_NS == p.row_NN && p.col_NS == p.col_NN ) {              if ( p.numSides == 1  ) {
317                 Finley_Assemble_PDE_Single2_3D(p,elements,S,F,A,B,C,D,X,Y);                 Finley_Assemble_PDE_Single2_3D(p,elements,S,F,A,B,C,D,X,Y);
318              } else if ( p.row_NS == p.col_NS &&  2*p.row_NS == p.row_NN && 2*p.col_NS == p.col_NN ) {              } else if ( p.numSides == 2 ) {
319                 if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {                 if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
320                    Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");                    Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
321                 } else {                 } else {
# Line 311  void Finley_Assemble_PDE(Finley_NodeFile Line 325  void Finley_Assemble_PDE(Finley_NodeFile
325                 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");                 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
326              }              }
327            } else if (p.numDim==2) {            } else if (p.numDim==2) {
328              if ((p.row_NS == p.col_NS) && (p.row_NS == p.row_NN) && (p.col_NS == p.col_NN )) {              if ( p.numSides == 1 ) {
329                 Finley_Assemble_PDE_Single2_2D(p,elements,S,F,A,B,C,D,X,Y);                 Finley_Assemble_PDE_Single2_2D(p,elements,S,F,A,B,C,D,X,Y);
330              } else if ( p.row_NS == p.col_NS &&  2*p.row_NS == p.row_NN && 2*p.col_NS == p.col_NN ) {              } else if ( p.numSides == 2 ) {
331                 if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {                 if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
332                    Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");                    Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
333                 } else {                 } else {
# Line 322  void Finley_Assemble_PDE(Finley_NodeFile Line 336  void Finley_Assemble_PDE(Finley_NodeFile
336              } else {              } else {
337                 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");                 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
338              }              }
339            } else if (p.numDim==2) {            } else if (p.numDim==1) {
340              if (p.row_NS == p.col_NS && p.row_NS == p.row_NN && p.col_NS == p.col_NN ) {              if ( p.numSides == 1 ) {
341                 Finley_Assemble_PDE_Single2_1D(p,elements,S,F,A,B,C,D,X,Y);                 Finley_Assemble_PDE_Single2_1D(p,elements,S,F,A,B,C,D,X,Y);
342              } else if ( p.row_NS == p.col_NS &&  2*p.row_NS == p.row_NN && 2*p.col_NS == p.col_NN ) {              } else if ( p.numSides == 2  ) {
343                 if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {                 if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
344                    Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");                    Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
345                 } else {                 } else {
# Line 341  void Finley_Assemble_PDE(Finley_NodeFile Line 355  void Finley_Assemble_PDE(Finley_NodeFile
355       } else {       } else {
356            Finley_setError(VALUE_ERROR,"Finley_Assemble_PDE requires number of equations == number of solutions  .");            Finley_setError(VALUE_ERROR,"Finley_Assemble_PDE requires number of equations == number of solutions  .");
357       }       }
      #ifdef Finley_TRACE  
      printf("timing: assemblage PDE: %.4e sec\n",Finley_timer()-time0);  
      #endif  
358    }    }
359      blocktimer_increment("Finley_Assemble_PDE()", blocktimer_start);
360  }  }
 /*  
  * $Log$  
  * Revision 1.8  2005/09/15 03:44:21  jgs  
  * Merge of development branch dev-02 back to main trunk on 2005-09-15  
  *  
  * Revision 1.7  2005/09/01 03:31:35  jgs  
  * Merge of development branch dev-02 back to main trunk on 2005-09-01  
  *  
  * Revision 1.6  2005/08/12 01:45:42  jgs  
  * erge of development branch dev-02 back to main trunk on 2005-08-12  
  *  
  * Revision 1.5.2.3  2005/09/07 06:26:17  gross  
  * the solver from finley are put into the standalone package paso now  
  *  
  * Revision 1.5.2.2  2005/08/24 02:02:18  gross  
  * timing output switched off. solver output can be swiched through getSolution(verbose=True) now.  
  *  
  * Revision 1.5.2.1  2005/08/03 08:54:27  gross  
  * contact element assemblage was called with wrong element table pointer  
  *  
  * Revision 1.5  2005/07/08 04:07:46  jgs  
  * Merge of development branch back to main trunk on 2005-07-08  
  *  
  * Revision 1.4  2004/12/15 07:08:32  jgs  
  * *** empty log message ***  
  * Revision 1.1.1.1.2.2  2005/06/29 02:34:47  gross  
  * some changes towards 64 integers in finley  
  *  
  * Revision 1.1.1.1.2.1  2004/11/24 01:37:12  gross  
  * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now  
  *  
  *  
  *  
  */  

Legend:
Removed from v.971  
changed lines
  Added in v.2748

  ViewVC Help
Powered by ViewVC 1.1.26