/[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.c revision 552 by gross, Wed Feb 22 03:18:07 2006 UTC trunk/finley/src/ElementFile_jacobeans.c revision 776 by gross, Wed Jul 12 00:07:31 2006 UTC
# Line 1  Line 1 
1  /*  /*
2   ******************************************************************************   ************************************************************
3   *                                                                            *   *          Copyright 2006 by ACcESS MNRF                   *
4   *       COPYRIGHT  ACcESS 2003,2004,2005,2006 -  All Rights Reserved              *   *                                                          *
5   *                                                                            *   *              http://www.access.edu.au                    *
6   * This software is the property of ACcESS. No part of this code              *   *       Primary Business: Queensland, Australia            *
7   * may be copied in any form or by any means without the expressed written    *   *  Licensed under the Open Software License version 3.0    *
8   * consent of ACcESS.  Copying, use or modification of this software          *   *     http://www.opensource.org/licenses/osl-3.0.php       *
9   * by any unauthorised person is illegal unless that person has a software    *   *                                                          *
10   * license agreement with ACcESS.                                             *   ************************************************************
  *                                                                            *  
  ******************************************************************************  
11  */  */
12    
13    
# Line 22  Line 20 
20    
21  /**************************************************************/  /**************************************************************/
22    
23  #include "ElementFile.h"  #include "Assemble.h"
 #include "Util.h"  
24  #ifdef _OPENMP  #ifdef _OPENMP
25  #include <omp.h>  #include <omp.h>
26  #endif  #endif
# Line 31  Line 28 
28    
29  /**************************************************************/  /**************************************************************/
30    
31  double* Finley_ElementFile_borrowDSDV(Finley_ElementFile* self) {  Finley_ElementFile_Jacobeans* Finley_ElementFile_Jacobeans_alloc(void)
32  }  {
33  double* Finley_ElementFile_borrowDSLinearDV(Finley_ElementFile* self) {    Finley_ElementFile_Jacobeans* out=MEMALLOC(1,Finley_ElementFile_Jacobeans);
34  }    if (Finley_checkPtr(out)) {
35  double* Finley_ElementFile_borrowLocalVolume(Finley_ElementFile* self) {       return NULL;
36      } else {
37     if (! self->volume_is_valid) {       out->status=FINLEY_INITIAL_STATUS-1;
38        numDim_t local_numDim= ;       out->volume=NULL;
39        numDim_t numDim= ;       out->DSDX=NULL;
40        numDim_t numQuad= ;       return out;
       if (self->volume==NULL) self->volume=MEMALLOC(self->numElements,double);  
       if (self->DvDV==NULL) self->DvDV=MEMALLOC(self->numElements*local_numDim*numDim,double);  
       if (! (Finley_checkPtr(self->volume) || Finley_checkPtr(self->DvDV)) ) {  
            self->volume_is_valid=TRUE;  
   
          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);  
41    }    }
42    }
43    
44    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);  
   }  
45    
46    if (! numSamplesEqual(C,p.numQuad,elements->numElements) ) {  void Finley_ElementFile_Jacobeans_dealloc(Finley_ElementFile_Jacobeans* in)
47          sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient C don't match (%d,%d)",p.numQuad,elements->numElements);  {
48          Finley_setError(TYPE_ERROR,error_msg);    if (in!=NULL) {
49        if (in->volume!=NULL) MEMFREE(in->volume);
50        if (in->DSDX!=NULL) MEMFREE(in->DSDX);  
51        MEMFREE(in);
52    }    }
53    }
54    
55    if (! numSamplesEqual(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);  
         Finley_setError(TYPE_ERROR,error_msg);  
   }  
56    
   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);  
   }  
57    
58    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);  
   }  
59    
60    /*  check the numDimensions: */  Finley_ElementFile_Jacobeans* Finley_ElementFile_borrowJacobeans(Finley_ElementFile* self, Finley_NodeFile* nodes,
61                                                                 bool_t reducedShapefunction, bool_t reducedIntegrationOrder) {
62      Finley_ElementFile_Jacobeans *out = NULL;
63        
64    if (p.numEqu==1 && p.numComp==1) {    if (reducedShapefunction) {
65      if (!isEmpty(A)) {         if (reducedIntegrationOrder) {
66        numDimensions[0]=p.numDim;             out=self->jacobeans_reducedS_reducedQ;
67        numDimensions[1]=p.numDim;         } else {
68        if (!isDataPointShapeEqual(A,2,numDimensions)) {             out=self->jacobeans_reducedS;
           sprintf(error_msg,"Finley_Assemble_PDE: coefficient A: illegal shape, expected shape (%d,%d)",numDimensions[0],numDimensions[1]);  
           Finley_setError(TYPE_ERROR,error_msg);  
       }  
     }  
     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);  
69         }         }
70      }    } else {
71      if (!isEmpty(D)) {         if (reducedIntegrationOrder) {
72         if (!isDataPointShapeEqual(D,0,numDimensions)) {             out=self->jacobeans_reducedQ;
73            Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficient D, rank 0 expected.");         } else {
74               out=self->jacobeans;
75         }         }
76      }    }
77      if (!isEmpty(X)) {  
78         numDimensions[0]=p.numDim;    if (out->status < nodes->status) {
79         if (!isDataPointShapeEqual(X,1,numDimensions)) {       Finley_RefElement *shape, *test;
80            sprintf(error_msg,"Finley_Assemble_PDE: coefficient X, expected shape (%d,",numDimensions[0]);       if (reducedShapefunction) {
81            Finley_setError(TYPE_ERROR,error_msg);         if (reducedIntegrationOrder) {
82               shape=self->LinearReferenceElement;
83               test=self->LinearReferenceElement;
84           } else {
85               shape=self->LinearReferenceElement;
86               test=self->LinearReferenceElement;
87         }         }
88      }       } else {
89      if (!isEmpty(Y)) {         if (reducedIntegrationOrder) {
90         if (!isDataPointShapeEqual(Y,0,numDimensions)) {             shape=self->ReferenceElement;
91            Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficient Y, rank 0 expected.");             test=self->ReferenceElement;
92           } else {
93               shape=self->ReferenceElement;
94               test=self->ReferenceElement;
95         }         }
96      }       }
   } else {  
     if (!isEmpty(A)) {  
       numDimensions[0]=p.numEqu;  
       numDimensions[1]=p.numDim;  
       numDimensions[2]=p.numComp;  
       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);  
       }  
     }  
   }  
97    
98    if (Finley_noError()) {       if (out->volume==NULL) out->volume=MEMALLOC((self->numElements)*(shape->numQuadNodes),double);
99       time0=Finley_timer();       if (out->DSDX==NULL) out->DSDX=MEMALLOC((self->numElements)*(test->Type->numShapes)*(shape->Type->numDim)*(shape->numQuadNodes),double);
100       #pragma omp parallel private(index_col,index_row,EM_S,EM_F,V,dVdv,dSdV,Vol,dvdV,color,q) \       if (! (Finley_checkPtr(out->volume) || Finley_checkPtr(out->DSDX)) ) {
101              firstprivate(elements,nodes,S,F,A,B,C,D,X,Y)            if (nodes->numDim==1) {
102       {               if (shape->Type->numDim==0) {
103           EM_S=EM_F=V=dVdv=dSdV=Vol=dvdV=NULL;  
104           index_row=index_col=NULL;               } else if (shape->Type->numDim==1) {
105                      Assemble_jacobeans_1D(nodes->Coordinates,shape->numQuadNodes,shape->QuadWeights,shape->Type->numShapes,
106           /* allocate work arrays: */                                          self->numElements,self->Nodes,
107                                            shape->dSdv,test->Type->numShapes,test->dSdv,
108           EM_S=(double*) THREAD_MEMALLOC(p.NN_row*p.NN_col*p.numEqu*p.numComp,double);                                          out->DSDX,out->volume,self->Id);
109           EM_F=(double*) THREAD_MEMALLOC(p.NN_row*p.numEqu,double);               } else {
110           V=(double*) THREAD_MEMALLOC(p.NN*p.numDim,double);                    Finley_setError(SYSTEM_ERROR,"Finley_ElementFile_borrowJacobeans: local dimenion in a 1D domain has to be less or equal 1.");
111           dVdv=(double*) THREAD_MEMALLOC(p.numDim*p.numDim*p.numQuad,double);               }
112           dvdV=(double*) THREAD_MEMALLOC(p.numDim*p.numDim*p.numQuad,double);            } else if (nodes->numDim==2) {
113           dSdV=(double*) THREAD_MEMALLOC(p.NS_row*p.numQuad*p.numDim,double);               if (shape->Type->numDim==0) {
114           Vol=(double*) THREAD_MEMALLOC(p.numQuad,double);  
115           index_col=(index_t*) THREAD_MEMALLOC(p.NN_col,index_t);               } else if (shape->Type->numDim==1) {
116           index_row=(index_t*) THREAD_MEMALLOC(p.NN_row,index_t);                    Assemble_jacobeans_2D_M1D(nodes->Coordinates,shape->numQuadNodes,shape->QuadWeights,shape->Type->numShapes,
117                                                self->numElements,self->Nodes,
118           if (! (Finley_checkPtr(EM_S) || Finley_checkPtr(EM_F) || Finley_checkPtr(V) || Finley_checkPtr(index_col) ||                                              shape->dSdv,test->Type->numShapes,test->dSdv,
119                  Finley_checkPtr(index_row) || Finley_checkPtr(dVdv) || Finley_checkPtr(dSdV) || Finley_checkPtr(Vol) ))  {                                              out->DSDX,out->volume,self->Id);
120                 } else if (shape->Type->numDim==2) {
121             /*  open loop over all colors: */                    Assemble_jacobeans_2D(nodes->Coordinates,shape->numQuadNodes,shape->QuadWeights,shape->Type->numShapes,
122             for (color=elements->minColor;color<=elements->maxColor;color++) {                                          self->numElements,self->Nodes,
123                /*  open loop over all elements: */                                          shape->dSdv,test->Type->numShapes,test->dSdv,
124                #pragma omp for private(e) schedule(static)                                          out->DSDX,out->volume,self->Id);
125                for(e=0;e<elements->numElements;e++){               } else {
126                  if (elements->Color[e]==color) {                 Finley_setError(SYSTEM_ERROR,"Finley_ElementFile_borrowJacobeans: local dimenion in a 2D domain has to be less or equal 2.");
127                    for (q=0;q<p.NN_row;q++) index_row[q]=p.label_row[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];               }
128                    /* gather V-coordinates of nodes into V: */            } else if (nodes->numDim==3) {
129  //============================               if (shape->Type->numDim==0) {
130            Finley_Util_Gather_double(NN,&(self->Nodes[INDEX2(0,e,NN)]),numDim,nodes->Coordinates,V);  
131                 } else if (shape->Type->numDim==1) {
132  void Finley_LocalVolume_33(ElementFile* elem,NodeFile* nodes) {                    Assemble_jacobeans_3D_M1D(nodes->Coordinates,shape->numQuadNodes,shape->QuadWeights,shape->Type->numShapes,
133      char error_msg[LenErrorMsg_MAX];                                              self->numElements,self->Nodes,
134      double tol=TOL*(node->diameter)**(node->dim);                                              shape->dSdv,test->Type->numShapes,test->dSdv,
135      #pragma omp parallel                                              out->DSDX,out->volume,self->Id);
136      {               } else if (shape->Type->numDim==2) {
137         register double DVDv00,DVDv10,DVDv20,DVDv01,DVDv11,DVDv21,DVDv02,DVDv12,DVDv22,D;                    Assemble_jacobeans_3D_M2D(nodes->Coordinates,shape->numQuadNodes,shape->QuadWeights,shape->Type->numShapes,
138         register double x0_tmp, x1_tmp, x2_tmp, dSdv0_tmp, dSdv1_tmp, dSdv2_tmp;                                              self->numElements,self->Nodes,
139         double x0[NN],x1[NN],x2[NN];                                              shape->dSdv,test->Type->numShapes,test->dSdv,
140         index_t q,k;                                              out->DSDX,out->volume,self->Id);
141         #pragma omp for private(e) schedule(static)               } else if (shape->Type->numDim==3) {
142         for(e=0;e<elements->numElements;e++){                    Assemble_jacobeans_3D(nodes->Coordinates,shape->numQuadNodes,shape->QuadWeights,shape->Type->numShapes,
143                                            self->numElements,self->Nodes,
144            for (k=0;k<NN;++k) {                                          shape->dSdv,test->Type->numShapes,test->dSdv,
145               x0[k]=nodes->Coordinates[INDEX2(0,elem->Nodes[INDEX2(k,e,NN)],numDim)];                                          out->DSDX,out->volume,self->Id);
146               x1[k]=nodes->Coordinates[INDEX2(1,elem->Nodes[INDEX2(k,e,NN)],numDim)];               } else {
147               x2[k]=nodes->Coordinates[INDEX2(2,elem->Nodes[INDEX2(k,e,NN)],numDim)];                 Finley_setError(SYSTEM_ERROR,"Finley_ElementFile_borrowJacobeans: local dimenion in a 3D domain has to be less or equal 3.");
148                 }
149              } else {
150                Finley_setError(SYSTEM_ERROR,"Finley_ElementFile_borrowJacobeans: spatial dimension has to be less or equal 3.");
151            }            }
152         }
153         if (Finley_noError()) {
154             out->status = nodes->status;
155         } else {
156             out=NULL;
157         }
158    
159            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;  
       
             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;  
     
             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];  
             }  
         }  
   
 }  
   
   
       
   
160    
161                    /*  calculate dVdv(i,j,q)=V(i,k)*DSDv(k,j,q) */     return out;
           Finley_Util_SmallMatMult(numDim,numDim*numQuad,dVdv,NN,V,referenceElement->dSdv);  
                   /*  dvdV=invert(dVdv) inplace */  
           Finley_Util_InvertSmallMat(numQuad,numDim,dVdv,dvdV,Vol);  
           for (q=0;q<numQuad;q++) self->volume[q]=ABS(self.volume[q]*p.referenceElement->QuadWeights[q]);  
                   /*  calculate dSdV=DSDv*DvDV */  
           Finley_Util_SmallMatSetMult(p.numQuad,p.NS_row,p.numDim,dSdV,p.numDim,p.referenceElement_row->dSdv,dvdV);  
                   /*  scale volume: */  
 //============================  
       
                    /*   integration for the stiffness matrix: */  
                    /*   in order to optimze the number of operations the case of constants coefficience needs a bit more work */  
                    /*   to make use of some symmetry. */  
   
                      if (S!=NULL) {  
                        for (q=0;q<p.NN_row*p.NN_col*p.numEqu*p.numComp;q++) EM_S[q]=0;  
                        if (p.numEqu==1 && p.numComp==1) {  
                        Finley_Assemble_PDEMatrix_Single2(p.NS_row,p.numDim,p.numQuad,  
                                                                      p.referenceElement_row->S,dSdV,Vol,p.NN_row,EM_S,  
                                                                      getSampleData(A,e),isExpanded(A),  
                                                                      getSampleData(B,e),isExpanded(B),  
                                                                      getSampleData(C,e),isExpanded(C),  
                                                                      getSampleData(D,e),isExpanded(D));  
                        } else {  
                        Finley_Assemble_PDEMatrix_System2(p.NS_row,p.numDim,p.numQuad,p.numEqu,p.numComp,  
                                                                      p.referenceElement_row->S,dSdV,Vol,p.NN_row,EM_S,  
                                                                      getSampleData(A,e),isExpanded(A),  
                                                                      getSampleData(B,e),isExpanded(B),  
                                                                      getSampleData(C,e),isExpanded(C),  
                                                                      getSampleData(D,e),isExpanded(D));  
                        }  
                        for (q=0;q<p.NN_col;q++) index_col[q]=p.label_col[elements->Nodes[INDEX2(p.col_node[q],e,p.NN)]];  
                        Finley_Assemble_addToSystemMatrix(S,p.NN_row,index_row,p.numEqu,p.NN_col,index_col,p.numComp,EM_S);  
                      }  
                      if (!isEmpty(F)) {  
                        for (q=0;q<p.NN_row*p.numEqu;q++) EM_F[q]=0;  
                        if (p.numEqu==1) {  
                        Finley_Assemble_RHSMatrix_Single(p.NS_row,p.numDim,p.numQuad,  
                                                                  p.referenceElement_row->S,dSdV,Vol,p.NN_row,EM_F,  
                                                                  getSampleData(X,e),isExpanded(X),  
                                                                  getSampleData(Y,e),isExpanded(Y));  
                        } else {  
                        Finley_Assemble_RHSMatrix_System(p.NS_row,p.numDim,p.numQuad,  
                                                                  p.numEqu,p.referenceElement_row->S,dSdV,Vol,p.NN_row,EM_F,  
                                                                  getSampleData(X,e),isExpanded(X),  
                                                                  getSampleData(Y,e),isExpanded(Y));  
                        }  
                        /* add  */  
                        Finley_Util_AddScatter(p.NN_row,index_row,p.numEqu,EM_F,getSampleData(F,0));  
                     }  
                 }  
               }  
             }  
          }  
          /* clean up */  
          THREAD_MEMFREE(EM_S);  
          THREAD_MEMFREE(EM_F);  
          THREAD_MEMFREE(V);  
          THREAD_MEMFREE(dVdv);  
          THREAD_MEMFREE(dvdV);  
          THREAD_MEMFREE(dSdV);  
          THREAD_MEMFREE(Vol);  
          THREAD_MEMFREE(index_col);  
          THREAD_MEMFREE(index_row);  
      }  
      #ifdef Finley_TRACE  
      printf("timing: assemblage PDE: %.4e sec\n",Finley_timer()-time0);  
      #endif  
   }  
162  }  }
163  /*  /*
164   * $Log$   * $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  
  *  
  *  
165   *   *
166   */   */

Legend:
Removed from v.552  
changed lines
  Added in v.776

  ViewVC Help
Powered by ViewVC 1.1.26