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

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

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

trunk/finley/src/ElementFile_borrowLocalVolume.cccc revision 700 by gross, Thu Apr 6 00:13:40 2006 UTC branches/domexper/dudley/src/ElementFile_jacobeans.c revision 3157 by jfenwick, Mon Sep 6 05:10:18 2010 UTC
# Line 1  Line 1 
 /*  
  ************************************************************  
  *          Copyright 2006 by ACcESS MNRF                   *  
  *                                                          *  
  *              http://www.access.edu.au                    *  
  *       Primary Business: Queensland, Australia            *  
  *  Licensed under the Open Software License version 3.0    *  
  *     http://www.opensource.org/licenses/osl-3.0.php       *  
  *                                                          *  
  ************************************************************  
 */  
1    
2    /*******************************************************
3    *
4    * Copyright (c) 2003-2010 by University of Queensland
5    * Earth Systems Science Computational Center (ESSCC)
6    * http://www.uq.edu.au/esscc
7    *
8    * 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    
 /**************************************************************/  
   
 /**************************************************************/  
   
 /*  Author: gross@access.edu.au */  
 /*  Version: $Id$ */  
   
 /**************************************************************/  
14    
15  #include "ElementFile.h"  #include "ElementFile.h"
16  #include "Util.h"  #include "Assemble.h"
17  #ifdef _OPENMP  #ifdef _OPENMP
18  #include <omp.h>  #include <omp.h>
19  #endif  #endif
# Line 29  Line 21 
21    
22  /**************************************************************/  /**************************************************************/
23    
24  double* Finley_ElementFile_borrowDSDV(Finley_ElementFile* self) {  Dudley_ElementFile_Jacobeans* Dudley_ElementFile_Jacobeans_alloc(Dudley_ShapeFunction* BasisFunctions)
25  }  {
26  double* Finley_ElementFile_borrowDSLinearDV(Finley_ElementFile* self) {    Dudley_ElementFile_Jacobeans* out=MEMALLOC(1,Dudley_ElementFile_Jacobeans);
27  }    if (Dudley_checkPtr(out)) {
28  double* Finley_ElementFile_borrowLocalVolume(Finley_ElementFile* self) {       return NULL;
29      } else {
30     if (! self->volume_is_valid) {       out->status=DUDLEY_INITIAL_STATUS-1;
31        numDim_t local_numDim= ;       out->BasisFunctions=Dudley_ShapeFunction_reference(BasisFunctions);
32        numDim_t numDim= ;       out->numDim=0;
33        numDim_t numQuad= ;       out->numQuadTotal=0;
34        if (self->volume==NULL) self->volume=MEMALLOC(self->numElements,double);       out->numElements=0;
35        if (self->DvDV==NULL) self->DvDV=MEMALLOC(self->numElements*local_numDim*numDim,double);       out->volume=NULL;
36        if (! (Finley_checkPtr(self->volume) || Finley_checkPtr(self->DvDV)) ) {       out->DSDX=NULL;
37             self->volume_is_valid=TRUE;       return out;
   
          if (numDim==1 && local_numDim==1) {  
          } else if (numDim==2 && local_numDim==1) {  
   
          } else if (numDim==2 && local_numDim==2) {  
   
          } else if (numDim==3 && local_numDim==1) {  
   
          } else if (numDim==3 && local_numDim==2) {  
           
          } else if (numDim==3 && local_numDim==3) {  
   
          }  
       }  
    }  
   
    return self->volume;  
 }  
   
   char error_msg[LenErrorMsg_MAX];  
   double *EM_S=NULL,*EM_F=NULL,*V=NULL,*dVdv=NULL,*dSdV=NULL,*Vol=NULL,*dvdV=NULL;  
   double time0;  
   numDim_t numDimensions[ESCRIPT_MAX_DATA_RANK],e,q;  
   Assemble_Parameters p;  
   index_t *index_row=NULL,*index_col=NULL,color;  
   Finley_resetError();  
   
   if (nodes==NULL || elements==NULL) return;  
   if (S==NULL && isEmpty(F)) return;  
   
   /* set all parameters in p*/  
   Assemble_getAssembleParameters(nodes,elements,S,F,&p);  
   if (! Finley_noError()) return;  
   
   /*  this function assumes NS=NN */  
   if (p.NN!=p.NS) {  
     Finley_setError(SYSTEM_ERROR,"Finley_Assemble_PDE: for Finley_Assemble_PDE numNodes and numShapes have to be identical.");  
     return;  
   }  
   if (p.numDim!=p.numElementDim) {  
     Finley_setError(SYSTEM_ERROR,"Finley_Assemble_PDE: Finley_Assemble_PDE accepts volume elements only.");  
     return;  
   }  
   /*  get a functionspace */  
   type_t funcspace=UNKNOWN;  
   updateFunctionSpaceType(funcspace,A);  
   updateFunctionSpaceType(funcspace,B);  
   updateFunctionSpaceType(funcspace,C);  
   updateFunctionSpaceType(funcspace,D);  
   updateFunctionSpaceType(funcspace,X);  
   updateFunctionSpaceType(funcspace,Y);  
   if (funcspace==UNKNOWN) return; /* all  data are empty */  
   
   /* check if all function spaces are the same */  
   
   if (! functionSpaceTypeEqual(funcspace,A) ) {  
         Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient A");  
   }  
   if (! functionSpaceTypeEqual(funcspace,B) ) {  
         Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient B");  
   }  
   if (! functionSpaceTypeEqual(funcspace,C) ) {  
         Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient C");  
   }  
   if (! functionSpaceTypeEqual(funcspace,D) ) {  
         Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient D");  
   }  
   if (! functionSpaceTypeEqual(funcspace,X) ) {  
         Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient X");  
   }  
   if (! functionSpaceTypeEqual(funcspace,Y) ) {  
         Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient Y");  
   }  
   
   /* check if all function spaces are the same */  
   
   if (! numSamplesEqual(A,p.numQuad,elements->numElements) ) {  
         sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient A don't match (%d,%d)",p.numQuad,elements->numElements);  
         Finley_setError(TYPE_ERROR,error_msg);  
   }  
   
   if (! numSamplesEqual(B,p.numQuad,elements->numElements) ) {  
         sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient B don't match (%d,%d)",p.numQuad,elements->numElements);  
         Finley_setError(TYPE_ERROR,error_msg);  
38    }    }
39    }
40    
41    if (! numSamplesEqual(C,p.numQuad,elements->numElements) ) {  /**************************************************************/
         sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient C don't match (%d,%d)",p.numQuad,elements->numElements);  
         Finley_setError(TYPE_ERROR,error_msg);  
   }  
42    
43    if (! numSamplesEqual(D,p.numQuad,elements->numElements) ) {  void Dudley_ElementFile_Jacobeans_dealloc(Dudley_ElementFile_Jacobeans* in)
44          sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient D don't match (%d,%d)",p.numQuad,elements->numElements);  {
45          Finley_setError(TYPE_ERROR,error_msg);    if (in!=NULL) {
46        Dudley_ShapeFunction_dealloc(in->BasisFunctions);  
47        MEMFREE(in->volume);
48        MEMFREE(in->DSDX);
49        MEMFREE(in);
50    }    }
51    }
52    
53    if (! numSamplesEqual(X,p.numQuad,elements->numElements) ) {  /**************************************************************/
         sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient X don't match (%d,%d)",p.numQuad,elements->numElements);  
         Finley_setError(TYPE_ERROR,error_msg);  
   }  
54    
   if (! numSamplesEqual(Y,p.numQuad,elements->numElements) ) {  
         sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient Y don't match (%d,%d)",p.numQuad,elements->numElements);  
         Finley_setError(TYPE_ERROR,error_msg);  
   }  
55    
56    /*  check the numDimensions: */  Dudley_ElementFile_Jacobeans* Dudley_ElementFile_borrowJacobeans(Dudley_ElementFile* self, Dudley_NodeFile* nodes,
57                                                                     bool_t reducedShapefunction, bool_t reducedIntegrationOrder) {
58      Dudley_ElementFile_Jacobeans *out = NULL;
59      Dudley_ShapeFunction *shape=NULL, *basis;
60      Dudley_ReferenceElement*  refElement=NULL;
61      double *dBdv;
62        
63    if (p.numEqu==1 && p.numComp==1) {    dim_t numNodes=self->numNodes;
64      if (!isEmpty(A)) {    
65        numDimensions[0]=p.numDim;    if (reducedShapefunction) {
66        numDimensions[1]=p.numDim;         if (reducedIntegrationOrder) {
67        if (!isDataPointShapeEqual(A,2,numDimensions)) {             out=self->jacobeans_reducedQ;
68            sprintf(error_msg,"Finley_Assemble_PDE: coefficient A: illegal shape, expected shape (%d,%d)",numDimensions[0],numDimensions[1]);         } else {
69            Finley_setError(TYPE_ERROR,error_msg);             out=self->jacobeans;
       }  
     }  
     if (!isEmpty(B)) {  
        numDimensions[0]=p.numDim;  
        if (!isDataPointShapeEqual(B,1,numDimensions)) {  
           sprintf(error_msg,"Finley_Assemble_PDE: coefficient B: illegal shape (%d,)",numDimensions[0]);  
           Finley_setError(TYPE_ERROR,error_msg);  
        }  
     }  
     if (!isEmpty(C)) {  
        numDimensions[0]=p.numDim;  
        if (!isDataPointShapeEqual(C,1,numDimensions)) {  
           sprintf(error_msg,"Finley_Assemble_PDE: coefficient C, expected shape (%d,)",numDimensions[0]);  
           Finley_setError(TYPE_ERROR,error_msg);  
        }  
     }  
     if (!isEmpty(D)) {  
        if (!isDataPointShapeEqual(D,0,numDimensions)) {  
           Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficient D, rank 0 expected.");  
        }  
     }  
     if (!isEmpty(X)) {  
        numDimensions[0]=p.numDim;  
        if (!isDataPointShapeEqual(X,1,numDimensions)) {  
           sprintf(error_msg,"Finley_Assemble_PDE: coefficient X, expected shape (%d,",numDimensions[0]);  
           Finley_setError(TYPE_ERROR,error_msg);  
        }  
     }  
     if (!isEmpty(Y)) {  
        if (!isDataPointShapeEqual(Y,0,numDimensions)) {  
           Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficient Y, rank 0 expected.");  
70         }         }
     }  
71    } else {    } else {
72      if (!isEmpty(A)) {         if (reducedIntegrationOrder) {
73        numDimensions[0]=p.numEqu;             out=self->jacobeans_reducedQ;
74        numDimensions[1]=p.numDim;         } else {
75        numDimensions[2]=p.numComp;             out=self->jacobeans;
76        numDimensions[3]=p.numDim;         }
       if (!isDataPointShapeEqual(A,4,numDimensions)) {  
           sprintf(error_msg,"Finley_Assemble_PDE: coefficient A, expected shape (%d,%d,%d,%d)",numDimensions[0],numDimensions[1],numDimensions[2],numDimensions[3]);  
           Finley_setError(TYPE_ERROR,error_msg);  
       }  
     }  
     if (!isEmpty(B)) {  
       numDimensions[0]=p.numEqu;  
       numDimensions[1]=p.numDim;  
       numDimensions[2]=p.numComp;  
       if (!isDataPointShapeEqual(B,3,numDimensions)) {  
           sprintf(error_msg,"Finley_Assemble_PDE: coefficient B, expected shape (%d,%d,%d)",numDimensions[0],numDimensions[1],numDimensions[2]);  
           Finley_setError(TYPE_ERROR,error_msg);  
       }  
     }  
     if (!isEmpty(C)) {  
       numDimensions[0]=p.numEqu;  
       numDimensions[1]=p.numComp;  
       numDimensions[2]=p.numDim;  
       if (!isDataPointShapeEqual(C,3,numDimensions)) {  
           sprintf(error_msg,"Finley_Assemble_PDE: coefficient C, expected shape (%d,%d,%d)",numDimensions[0],numDimensions[1],numDimensions[2]);  
           Finley_setError(TYPE_ERROR,error_msg);  
       }  
     }  
     if (!isEmpty(D)) {  
       numDimensions[0]=p.numEqu;  
       numDimensions[1]=p.numComp;  
       if (!isDataPointShapeEqual(D,2,numDimensions)) {  
           sprintf(error_msg,"Finley_Assemble_PDE: coefficient D, expected shape (%d,%d)",numDimensions[0],numDimensions[1]);  
           Finley_setError(TYPE_ERROR,error_msg);  
       }  
     }  
     if (!isEmpty(X)) {  
       numDimensions[0]=p.numEqu;  
       numDimensions[1]=p.numDim;  
       if (!isDataPointShapeEqual(X,2,numDimensions)) {  
           sprintf(error_msg,"Finley_Assemble_PDE: coefficient X, expected shape (%d,%d)",numDimensions[0],numDimensions[1]);  
           Finley_setError(TYPE_ERROR,error_msg);  
       }  
     }  
     if (!isEmpty(Y)) {  
       numDimensions[0]=p.numEqu;  
       if (!isDataPointShapeEqual(Y,1,numDimensions)) {  
           sprintf(error_msg,"Finley_Assemble_PDE: coefficient Y, expected shape (%d,)",numDimensions[0]);  
           Finley_setError(TYPE_ERROR,error_msg);  
       }  
     }  
77    }    }
78      if (out->status < nodes->status) {
   if (Finley_noError()) {  
      time0=Finley_timer();  
      #pragma omp parallel private(index_col,index_row,EM_S,EM_F,V,dVdv,dSdV,Vol,dvdV,color,q) \  
             firstprivate(elements,nodes,S,F,A,B,C,D,X,Y)  
      {  
          EM_S=EM_F=V=dVdv=dSdV=Vol=dvdV=NULL;  
          index_row=index_col=NULL;  
   
          /* allocate work arrays: */  
   
          EM_S=(double*) THREAD_MEMALLOC(p.NN_row*p.NN_col*p.numEqu*p.numComp,double);  
          EM_F=(double*) THREAD_MEMALLOC(p.NN_row*p.numEqu,double);  
          V=(double*) THREAD_MEMALLOC(p.NN*p.numDim,double);  
          dVdv=(double*) THREAD_MEMALLOC(p.numDim*p.numDim*p.numQuad,double);  
          dvdV=(double*) THREAD_MEMALLOC(p.numDim*p.numDim*p.numQuad,double);  
          dSdV=(double*) THREAD_MEMALLOC(p.NS_row*p.numQuad*p.numDim,double);  
          Vol=(double*) THREAD_MEMALLOC(p.numQuad,double);  
          index_col=(index_t*) THREAD_MEMALLOC(p.NN_col,index_t);  
          index_row=(index_t*) THREAD_MEMALLOC(p.NN_row,index_t);  
   
          if (! (Finley_checkPtr(EM_S) || Finley_checkPtr(EM_F) || Finley_checkPtr(V) || Finley_checkPtr(index_col) ||  
                 Finley_checkPtr(index_row) || Finley_checkPtr(dVdv) || Finley_checkPtr(dSdV) || Finley_checkPtr(Vol) ))  {  
   
            /*  open loop over all colors: */  
            for (color=elements->minColor;color<=elements->maxColor;color++) {  
               /*  open loop over all elements: */  
               #pragma omp for private(e) schedule(static)  
               for(e=0;e<elements->numElements;e++){  
                 if (elements->Color[e]==color) {  
                   for (q=0;q<p.NN_row;q++) index_row[q]=p.label_row[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];  
                   /* gather V-coordinates of nodes into V: */  
 //============================  
           Finley_Util_Gather_double(NN,&(self->Nodes[INDEX2(0,e,NN)]),numDim,nodes->Coordinates,V);  
   
 void Finley_LocalVolume_33(ElementFile* elem,NodeFile* nodes) {  
     char error_msg[LenErrorMsg_MAX];  
     double tol=TOL*(node->diameter)**(node->dim);  
     #pragma omp parallel  
     {  
        register double DVDv00,DVDv10,DVDv20,DVDv01,DVDv11,DVDv21,DVDv02,DVDv12,DVDv22,D;  
        register double x0_tmp, x1_tmp, x2_tmp, dSdv0_tmp, dSdv1_tmp, dSdv2_tmp;  
        double x0[NN],x1[NN],x2[NN];  
        index_t q,k;  
        #pragma omp for private(e) schedule(static)  
        for(e=0;e<elements->numElements;e++){  
   
           for (k=0;k<NN;++k) {  
              x0[k]=nodes->Coordinates[INDEX2(0,elem->Nodes[INDEX2(k,e,NN)],numDim)];  
              x1[k]=nodes->Coordinates[INDEX2(1,elem->Nodes[INDEX2(k,e,NN)],numDim)];  
              x2[k]=nodes->Coordinates[INDEX2(2,elem->Nodes[INDEX2(k,e,NN)],numDim)];  
           }  
   
           for (q=0;q<numQuad;++q) {  
             x0_tmp=x0[0];  
             x1_tmp=x1[0];  
             x2_tmp=x2[0];  
     
             dSdv0_tmp=elem->referenceElement->dSdv[INDEX3(0,0,q,NN,local_numDim)];  
             dSdv1_tmp=elem->referenceElement->dSdv[INDEX3(0,1,q,NN,local_numDim)];  
             dSdv2_tmp=elem->referenceElement->dSdv[INDEX3(0,2,q,NN,local_numDim)];  
   
             DVDv00=x0_tmp*dSdv0_tmp;  
             DVDv10=x1_tmp*dSdv0_tmp;  
             DVDv20=x2_tmp*dSdv0_tmp;  
79            
80              DVDv01=x0_tmp*dSdv1_tmp;       basis=out->BasisFunctions;
81              DVDv11=x1_tmp*dSdv1_tmp;       shape=Dudley_ReferenceElementSet_borrowParametrization(self->referenceElementSet, reducedIntegrationOrder);
82              DVDv21=x2_tmp*dSdv1_tmp;       refElement= Dudley_ReferenceElementSet_borrowReferenceElement(self->referenceElementSet, reducedIntegrationOrder);
83      
84              DVDv02=x0_tmp*dSdv2_tmp;       out->numDim=nodes->numDim;
85              DVDv12=x1_tmp*dSdv2_tmp;       out->numQuadTotal=shape->numQuadNodes;
86              DVDv22=x2_tmp*dSdv2_tmp;       out->numShapesTotal=basis->Type->numShapes;
87           out->numElements=self->numElements;
             for (k=1;k<NN;++k) {  
                 x0_tmp=x0[k];  
                 x1_tmp=x1[k];  
                 x2_tmp=x2[k];  
         
                 dSdv0_tmp=elem->referenceElement->dSdv[INDEX3(k,0,q,NN,local_numDim)];  
                 dSdv1_tmp=elem->referenceElement->dSdv[INDEX3(k,1,q,NN,local_numDim)];  
                 dSdv2_tmp=elem->referenceElement->dSdv[INDEX3(k,2,q,NN,local_numDim)];  
         
                 DVDv00+=x0_tmp*dSdv0_tmp;  
                 DVDv10+=x1_tmp*dSdv0_tmp;  
                 DVDv20+=x2_tmp*dSdv0_tmp;  
           
                 DVDv01+=x0_tmp*dSdv1_tmp;  
                 DVDv11+=x1_tmp*dSdv1_tmp;  
                 DVDv21+=x2_tmp*dSdv1_tmp;  
         
                 DVDv02+=x0_tmp*dSdv2_tmp;  
                 DVDv12+=x1_tmp*dSdv2_tmp;  
                 DVDv22+=x2_tmp*dSdv2_tmp;  
             }  
             D = DVDv00*(DVDv11*DVDv22-DVDv12*DVDv21)+  
                 DVDv01*(DVDv20*DVDv12-DVDv10*DVDv22)+  
                 DVDv02*(DVDv10*DVDv21-DVDv20*DVDv11);  
             if (ABS(D)<tol){  
                 sprintf(error_msg,"Finley_LocalVolume_33: volume element %d of type %s  ",e,elem->referenceElement->name,ABS(D),tol);  
                 Finley_setError(ZERO_DIVISION_ERROR,error_msg);  
             } else {  
                 D=1./D;  
                 elem->DvDV[INDEX3(0,0,q,local_numDim,numDim)]=(DVDv11*DVDv22-DVDv12*DVDv21)*D;  
                 elem->DvDV[INDEX3(1,0,q,local_numDim,numDim)]=(DVDv20*DVDv12-DVDv10*DVDv22)*D;  
                 elem->DvDV[INDEX3(2,0,q,local_numDim,numDim)]=(DVDv10*DVDv21-DVDv20*DVDv11)*D;  
                 elem->DvDV[INDEX3(0,1,q,local_numDim,numDim)]=(DVDv02*DVDv21-DVDv01*DVDv22)*D;  
                 elem->DvDV[INDEX3(1,1,q,local_numDim,numDim)]=(DVDv00*DVDv22-DVDv20*DVDv02)*D;  
                 elem->DvDV[INDEX3(2,1,q,local_numDim,numDim)]=(DVDv01*DVDv20-DVDv00*DVDv21)*D;  
                 elem->DvDV[INDEX3(0,2,q,local_numDim,numDim)]=(DVDv01*DVDv12-DVDv02*DVDv11)*D;  
                 elem->DvDV[INDEX3(1,2,q,local_numDim,numDim)]=(DVDv02*DVDv10-DVDv00*DVDv12)*D;  
                 elem->DvDV[INDEX3(2,2,q,local_numDim,numDim)]=(DVDv00*DVDv11-DVDv01*DVDv10)*D;  
                 elem->volume[q]=ABS(D)*referenceElement->QuadWeights[q];  
             }  
         }  
   
 }  
   
   
88            
89         if (reducedShapefunction) {
90            dBdv=basis->dSdv;
91         } else {
92            dBdv=refElement->DBasisFunctionDv;
93         }
94        
95         if (out->numQuadTotal != basis->numQuadNodes) {
96            Dudley_setError(SYSTEM_ERROR,"Dudley_ElementFile_borrowJacobeans: Incorrect total number of quadrature points.");
97            return NULL;
98         }
99         if (refElement->numNodes> numNodes) {
100               Dudley_setError(SYSTEM_ERROR,"Dudley_ElementFile_borrowJacobeans: Too many nodes expected.");
101               return NULL;
102         }
103    
104         if (out->volume==NULL) out->volume=MEMALLOC((out->numElements)*(out->numQuadTotal),double);
105                    /*  calculate dVdv(i,j,q)=V(i,k)*DSDv(k,j,q) */       if (out->DSDX==NULL) out->DSDX=MEMALLOC((out->numElements)
106            Finley_Util_SmallMatMult(numDim,numDim*numQuad,dVdv,NN,V,referenceElement->dSdv);                                              *(out->numShapesTotal)
107                    /*  dvdV=invert(dVdv) inplace */                                              *(out->numDim)
108            Finley_Util_InvertSmallMat(numQuad,numDim,dVdv,dvdV,Vol);                                              *(out->numQuadTotal),double);
109            for (q=0;q<numQuad;q++) self->volume[q]=ABS(self.volume[q]*p.referenceElement->QuadWeights[q]);       if (! (Dudley_checkPtr(out->volume) || Dudley_checkPtr(out->DSDX)) ) {
110                    /*  calculate dSdV=DSDv*DvDV */            /*========================== dim = 1 ============================================== */
111            Finley_Util_SmallMatSetMult(p.numQuad,p.NS_row,p.numDim,dSdV,p.numDim,p.referenceElement_row->dSdv,dvdV);            if (out->numDim==1) {
112                    /*  scale volume: */               if (refElement->numLocalDim==0) {
113  //============================            Dudley_setError(SYSTEM_ERROR,"Dudley_ElementFile_borrowJacobeans: 1D does not support local dimension 0.");
114                     } else if (refElement->numLocalDim==1) {
115                     /*   integration for the stiffness matrix: */                         Assemble_jacobeans_1D(nodes->Coordinates,out->numQuadTotal,shape->QuadWeights,
116                     /*   in order to optimze the number of operations the case of constants coefficience needs a bit more work */                                              shape->Type->numShapes,self->numElements,numNodes,self->Nodes,
117                     /*   to make use of some symmetry. */                                              shape->dSdv,basis->Type->numShapes,dBdv, out->DSDX,out->volume,self->Id);
118    
119                       if (S!=NULL) {               } else {
120                         for (q=0;q<p.NN_row*p.NN_col*p.numEqu*p.numComp;q++) EM_S[q]=0;                    Dudley_setError(SYSTEM_ERROR,"Dudley_ElementFile_borrowJacobeans: local dimenion in a 1D domain has to be 0 or 1.");
121                         if (p.numEqu==1 && p.numComp==1) {               }
122                         Finley_Assemble_PDEMatrix_Single2(p.NS_row,p.numDim,p.numQuad,            /*========================== dim = 2 ============================================== */
123                                                                       p.referenceElement_row->S,dSdV,Vol,p.NN_row,EM_S,            } else if (out->numDim==2) {
124                                                                       getSampleData(A,e),isExpanded(A),               if (refElement->numLocalDim==0) {
125                                                                       getSampleData(B,e),isExpanded(B),                   Dudley_setError(SYSTEM_ERROR,"Dudley_ElementFile_borrowJacobeans: 2D does not support local dimension 0.");
126                                                                       getSampleData(C,e),isExpanded(C),               } else if (refElement->numLocalDim==1) {
127                                                                       getSampleData(D,e),isExpanded(D));                    if (out->BasisFunctions->Type->numDim==2) {
128                         } else {                          Assemble_jacobeans_2D_M1D_E2D(nodes->Coordinates,out->numQuadTotal,shape->QuadWeights,
129                         Finley_Assemble_PDEMatrix_System2(p.NS_row,p.numDim,p.numQuad,p.numEqu,p.numComp,                                                        shape->Type->numShapes,self->numElements,numNodes,self->Nodes,
130                                                                       p.referenceElement_row->S,dSdV,Vol,p.NN_row,EM_S,                                                        shape->dSdv,basis->Type->numShapes,dBdv,
131                                                                       getSampleData(A,e),isExpanded(A),                                                        out->DSDX,out->volume,self->Id);
132                                                                       getSampleData(B,e),isExpanded(B),                    }  else if (out->BasisFunctions->Type->numDim==1) {
133                                                                       getSampleData(C,e),isExpanded(C),                          Assemble_jacobeans_2D_M1D_E1D(nodes->Coordinates,out->numQuadTotal,shape->QuadWeights,
134                                                                       getSampleData(D,e),isExpanded(D));                                                        shape->Type->numShapes,self->numElements,numNodes,self->Nodes,
135                         }                                                        shape->dSdv,basis->Type->numShapes,dBdv,
136                         for (q=0;q<p.NN_col;q++) index_col[q]=p.label_col[elements->Nodes[INDEX2(p.col_node[q],e,p.NN)]];                                                        out->DSDX,out->volume,self->Id);
137                         Finley_Assemble_addToSystemMatrix(S,p.NN_row,index_row,p.numEqu,p.NN_col,index_col,p.numComp,EM_S);                    } else {
138                       }                      Dudley_setError(SYSTEM_ERROR,"Dudley_ElementFile_borrowJacobeans: element dimension for local dimenion 1 in a 2D domain has to be 1 or 2.");
139                       if (!isEmpty(F)) {                    }
140                         for (q=0;q<p.NN_row*p.numEqu;q++) EM_F[q]=0;               } else if (refElement->numLocalDim==2) {
141                         if (p.numEqu==1) {                       Assemble_jacobeans_2D(nodes->Coordinates,out->numQuadTotal,shape->QuadWeights,
142                         Finley_Assemble_RHSMatrix_Single(p.NS_row,p.numDim,p.numQuad,                                             shape->Type->numShapes,self->numElements,numNodes,self->Nodes,
143                                                                   p.referenceElement_row->S,dSdV,Vol,p.NN_row,EM_F,                                             shape->dSdv,basis->Type->numShapes,dBdv,
144                                                                   getSampleData(X,e),isExpanded(X),                                             out->DSDX,out->volume,self->Id);
145                                                                   getSampleData(Y,e),isExpanded(Y));               } else {
146                         } else {                 Dudley_setError(SYSTEM_ERROR,"Dudley_ElementFile_borrowJacobeans: local dimenion in a 2D domain has to be  1 or 2.");
147                         Finley_Assemble_RHSMatrix_System(p.NS_row,p.numDim,p.numQuad,               }
148                                                                   p.numEqu,p.referenceElement_row->S,dSdV,Vol,p.NN_row,EM_F,            /*========================== dim = 3 ============================================== */
149                                                                   getSampleData(X,e),isExpanded(X),            } else if (out->numDim==3) {
150                                                                   getSampleData(Y,e),isExpanded(Y));               if (refElement->numLocalDim==0) {
151                         }            Dudley_setError(SYSTEM_ERROR,"Dudley_ElementFile_borrowJacobeans: 3D does not support local dimension 0.");
152                         /* add  */               } else if (refElement->numLocalDim==2) {
153                         Finley_Util_AddScatter(p.NN_row,index_row,p.numEqu,EM_F,getSampleData(F,0));                    if (out->BasisFunctions->Type->numDim==3) {
154                      }                          Assemble_jacobeans_3D_M2D_E3D(nodes->Coordinates,out->numQuadTotal,shape->QuadWeights,
155                  }                                                        shape->Type->numShapes,self->numElements,numNodes,self->Nodes,
156                }                                                        shape->dSdv,basis->Type->numShapes,dBdv,
157              }                                                        out->DSDX,out->volume,self->Id);
158           }                    }  else if (out->BasisFunctions->Type->numDim==2) {
159           /* clean up */                          Assemble_jacobeans_3D_M2D_E2D(nodes->Coordinates,out->numQuadTotal,shape->QuadWeights,
160           THREAD_MEMFREE(EM_S);                                                        shape->Type->numShapes,self->numElements,numNodes,self->Nodes,
161           THREAD_MEMFREE(EM_F);                                                        shape->dSdv,basis->Type->numShapes,dBdv,
162           THREAD_MEMFREE(V);                                                        out->DSDX,out->volume,self->Id);
163           THREAD_MEMFREE(dVdv);                    } else {
164           THREAD_MEMFREE(dvdV);                      Dudley_setError(SYSTEM_ERROR,"Dudley_ElementFile_borrowJacobeans: element dimension for local dimenion 2 in a 3D domain has to be 3 or 2.");
165           THREAD_MEMFREE(dSdV);                    }
166           THREAD_MEMFREE(Vol);               } else if (refElement->numLocalDim==3) {
167           THREAD_MEMFREE(index_col);                       Assemble_jacobeans_3D(nodes->Coordinates,out->numQuadTotal,shape->QuadWeights,
168           THREAD_MEMFREE(index_row);                                             shape->Type->numShapes,self->numElements,numNodes,self->Nodes,
169                                               shape->dSdv,basis->Type->numShapes,dBdv,
170                                               out->DSDX,out->volume,self->Id);
171                 } else {
172                   Dudley_setError(SYSTEM_ERROR,"Dudley_ElementFile_borrowJacobeans: local dimenion in a 3D domain has to be 2 or 3.");
173                 }
174              } else {
175                Dudley_setError(SYSTEM_ERROR,"Dudley_ElementFile_borrowJacobeans: spatial dimension has to be 1, 2 or 3.");
176              }
177       }       }
178       #ifdef Finley_TRACE       if (Dudley_noError()) {
179       printf("timing: assemblage PDE: %.4e sec\n",Finley_timer()-time0);           out->status = nodes->status;
180       #endif       } else {
181    }           out=NULL;
182         }
183    
184       }
185    
186       return out;
187  }  }
 /*  
  * $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.700  
changed lines
  Added in v.3157

  ViewVC Help
Powered by ViewVC 1.1.26