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

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

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

revision 964 by gross, Tue Feb 13 05:10:26 2007 UTC revision 1028 by gross, Wed Mar 14 00:15:24 2007 UTC
# Line 120  void Finley_Quad_getNodesTet(int numQuad Line 120  void Finley_Quad_getNodesTet(int numQuad
120  void Finley_Quad_getNodesRec(int numQuadNodes,double* quadNodes,double* quadWeights) {  void Finley_Quad_getNodesRec(int numQuadNodes,double* quadNodes,double* quadWeights) {
121    char error_msg[LenErrorMsg_MAX];    char error_msg[LenErrorMsg_MAX];
122    int numQuadNodes1d,i,j,l;    int numQuadNodes1d,i,j,l;
123    double quadNodes1d[numQuadNodes],quadWeights1d[numQuadNodes];    double *quadNodes1d=NULL,*quadWeights1d=NULL;
124      bool_t set=FALSE;
125    #define DIM 2    #define DIM 2
126        
127    /*  find numQuadNodes1d with numQuadNodes1d**2==numQuadNodes: */    quadNodes1d=TMPMEMALLOC(numQuadNodes,double);
128        quadWeights1d=TMPMEMALLOC(numQuadNodes,double);
129    for (numQuadNodes1d=1;numQuadNodes1d<=MAX_numQuadNodesLine;numQuadNodes1d++) {    if (! ( Finley_checkPtr(quadNodes1d) || Finley_checkPtr(quadWeights1d) ) ) {
130      if (numQuadNodes1d*numQuadNodes1d==numQuadNodes) {       /*  find numQuadNodes1d with numQuadNodes1d**2==numQuadNodes: */
131              
132        /*  get 1D scheme: */       for (numQuadNodes1d=1;numQuadNodes1d<=MAX_numQuadNodesLine;numQuadNodes1d++) {
133           if (numQuadNodes1d*numQuadNodes1d==numQuadNodes) {
134                
135        Finley_Quad_getNodesLine(numQuadNodes1d,quadNodes1d,quadWeights1d);           /*  get 1D scheme: */
136            
137             Finley_Quad_getNodesLine(numQuadNodes1d,quadNodes1d,quadWeights1d);
138                
139        /*  make 2D scheme: */           /*  make 2D scheme: */
140                
141        if (Finley_noError()) {           if (Finley_noError()) {
142          l=0;             l=0;
143          for (i=0;i<numQuadNodes1d;i++) {             for (i=0;i<numQuadNodes1d;i++) {
144            for (j=0;j<numQuadNodes1d;j++) {               for (j=0;j<numQuadNodes1d;j++) {
145              QUADNODES(0,l)=quadNodes1d[i];                 QUADNODES(0,l)=quadNodes1d[i];
146              QUADNODES(1,l)=quadNodes1d[j];                 QUADNODES(1,l)=quadNodes1d[j];
147              QUADWEIGHTS(l)=quadWeights1d[i]*quadWeights1d[j];                 QUADWEIGHTS(l)=quadWeights1d[i]*quadWeights1d[j];
148              l++;                 l++;
149            }               }
150          }             }
151        }             set=TRUE;
152        return;             break;
153      }           }
154    }         }
155    sprintf(error_msg,"Finley_Quad_getNodesRec: Illegal number of quadrature nodes %d on hexahedron.",numQuadNodes);       }
156    Finley_setError(VALUE_ERROR,error_msg);       if (!set) {
157    #undef DIM           sprintf(error_msg,"Finley_Quad_getNodesRec: Illegal number of quadrature nodes %d on hexahedron.",numQuadNodes);
158             Finley_setError(VALUE_ERROR,error_msg);
159         }
160         TMPMEMFREE(quadNodes1d);
161         TMPMEMFREE(quadWeights1d);
162       }
163       #undef DIM
164  }  }
165    
166  /**************************************************************/  /**************************************************************/
# Line 161  void Finley_Quad_getNodesRec(int numQuad Line 171  void Finley_Quad_getNodesRec(int numQuad
171  void Finley_Quad_getNodesHex(int numQuadNodes,double* quadNodes,double* quadWeights) {  void Finley_Quad_getNodesHex(int numQuadNodes,double* quadNodes,double* quadWeights) {
172    char error_msg[LenErrorMsg_MAX];    char error_msg[LenErrorMsg_MAX];
173    int numQuadNodes1d,i,j,k,l;    int numQuadNodes1d,i,j,k,l;
174    double quadNodes1d[numQuadNodes],quadWeights1d[numQuadNodes];    double *quadNodes1d=NULL,*quadWeights1d=NULL;
175      bool_t set;
176    #define DIM 3    #define DIM 3
177        
178    /*  find numQuadNodes1d with numQuadNodes1d**3==numQuadNodes: */    /*  find numQuadNodes1d with numQuadNodes1d**3==numQuadNodes: */
179        
180    for (numQuadNodes1d=1;numQuadNodes1d<=MAX_numQuadNodesLine;numQuadNodes1d++) {    quadNodes1d=TMPMEMALLOC(numQuadNodes,double);
181      if (numQuadNodes1d*numQuadNodes1d*numQuadNodes1d==numQuadNodes) {    quadWeights1d=TMPMEMALLOC(numQuadNodes,double);
182      if (! ( Finley_checkPtr(quadNodes1d) || Finley_checkPtr(quadWeights1d) ) ) {
183         for (numQuadNodes1d=1;numQuadNodes1d<=MAX_numQuadNodesLine;numQuadNodes1d++) {
184           if (numQuadNodes1d*numQuadNodes1d*numQuadNodes1d==numQuadNodes) {
185                
186        /*  get 1D scheme: */           /*  get 1D scheme: */
187                
188        Finley_Quad_getNodesLine(numQuadNodes1d,quadNodes1d,quadWeights1d);           Finley_Quad_getNodesLine(numQuadNodes1d,quadNodes1d,quadWeights1d);
189                
190        /*  make 3D scheme: */           /*  make 3D scheme: */
191                
192        if (Finley_noError()) {           if (Finley_noError()) {
193          l=0;             l=0;
194          for (i=0;i<numQuadNodes1d;i++) {             for (i=0;i<numQuadNodes1d;i++) {
195            for (j=0;j<numQuadNodes1d;j++) {               for (j=0;j<numQuadNodes1d;j++) {
196              for (k=0;k<numQuadNodes1d;k++) {                 for (k=0;k<numQuadNodes1d;k++) {
197                QUADNODES(0,l)=quadNodes1d[i];                   QUADNODES(0,l)=quadNodes1d[i];
198                QUADNODES(1,l)=quadNodes1d[j];                   QUADNODES(1,l)=quadNodes1d[j];
199                QUADNODES(2,l)=quadNodes1d[k];                   QUADNODES(2,l)=quadNodes1d[k];
200                QUADWEIGHTS(l)=quadWeights1d[i]*quadWeights1d[j]*quadWeights1d[k];                   QUADWEIGHTS(l)=quadWeights1d[i]*quadWeights1d[j]*quadWeights1d[k];
201                l++;                   l++;
202              }                 }
203            }               }
204          }             }
205        }             set=TRUE;
206                     break;
207        return;           }
208      }         }
209         }
210         if (!set) {
211              sprintf(error_msg,"Finley_Quad_getNodesHex: Illegal number of quadrature nodes %d on hexahedron.",numQuadNodes);
212              Finley_setError(VALUE_ERROR,error_msg);
213         }
214         TMPMEMFREE(quadNodes1d);
215         TMPMEMFREE(quadWeights1d);
216    }    }
   sprintf(error_msg,"Finley_Quad_getNodesHex: Illegal number of quadrature nodes %d on hexahedron.",numQuadNodes);  
   Finley_setError(VALUE_ERROR,error_msg);  
217    #undef DIM    #undef DIM
218  }  }
219    
# Line 385  void Finley_Quad_getNodesPointOnFace(int Line 404  void Finley_Quad_getNodesPointOnFace(int
404    
405  void Finley_Quad_makeNodesOnFace(int dim, int numQuadNodes,double* quadNodes,double* quadWeights, Finley_Quad_getNodes getFaceNodes) {  void Finley_Quad_makeNodesOnFace(int dim, int numQuadNodes,double* quadNodes,double* quadWeights, Finley_Quad_getNodes getFaceNodes) {
406      int q,i;      int q,i;
407      double quadNodesOnFace[numQuadNodes*(dim-1)];      double *quadNodesOnFace=NULL;
   
408      #define DIM dim      #define DIM dim
409      getFaceNodes(numQuadNodes,quadNodesOnFace,quadWeights);      quadNodesOnFace=TMPMEMALLOC(numQuadNodes*(dim-1),double);
410      if (! Finley_noError()) return;  
411            if (! Finley_checkPtr(quadNodesOnFace) ) {
412      for (q=0;q<numQuadNodes;q++) {         getFaceNodes(numQuadNodes,quadNodesOnFace,quadWeights);
413         for (i=0;i<dim-1;i++) QUADNODES(i,q)=quadNodesOnFace[INDEX2(i,q,dim-1)];         if (Finley_noError()) {
414         QUADNODES(dim-1,q)=0;            for (q=0;q<numQuadNodes;q++) {
415                 for (i=0;i<dim-1;i++) QUADNODES(i,q)=quadNodesOnFace[INDEX2(i,q,dim-1)];
416                 QUADNODES(dim-1,q)=0;
417              }
418           }
419           TMPMEMFREE(quadNodesOnFace);
420      }      }
421      #undef DIM      #undef DIM
422  }  }

Legend:
Removed from v.964  
changed lines
  Added in v.1028

  ViewVC Help
Powered by ViewVC 1.1.26