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

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

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

trunk/finley/src/Assemble_PDE.c revision 751 by bcumming, Mon Jun 26 01:46:34 2006 UTC trunk/finley/src/Assemble_jacobeans.c revision 776 by gross, Wed Jul 12 00:07:31 2006 UTC
# Line 13  Line 13 
13    
14  /**************************************************************/  /**************************************************************/
15    
 /*    assembles the system of numEq PDEs into the stiffness matrix S and right hand side F */  
   
 /*     -div(A*grad u)-div(B*u)+C*grad u + D*u= -div X + Y */  
   
 /*      -(A_{k,i,m,j} u_m,j)_i-(B_{k,i,m} u_m)_i+C_{k,m,j} u_m,j-D_{k,m} u_m = -(X_{k,i})_i + Y_k */  
   
 /*    u has numComp components. */  
   
 /*    Shape of the coefficients: */  
   
 /*      A = numEqu x numDim x numComp x numDim */  
 /*      B = numDim x numEqu x numComp  */  
 /*      C = numEqu x numDim x numComp  */  
 /*      D = numEqu x numComp  */  
 /*      X = numEqu x numDim   */  
 /*      Y = numEqu */  
   
 /*    The coefficients A,B,C,D,X and Y have to be defined on the integartion points or not present (=NULL). */  
   
 /*    S and F have to be initialized before the routine is called. S or F can be NULL. In this case the left or */  
 /*    the right hand side of the PDE is not processed.  */  
   
 /*    The routine does not consider any boundary conditions. */  
   
 /**************************************************************/  
   
16  /*  Author: gross@access.edu.au */  /*  Author: gross@access.edu.au */
17  /*  Version: $Id$ */  /*  Version: $Id$ */
18    
# Line 50  Line 24 
24  #include <omp.h>  #include <omp.h>
25  #endif  #endif
26    
   
27  /**************************************************************/  /**************************************************************/
28    /*                                                            */
29    /*  Jacobean 1D                                               */
30    /*                                                            */
31    void Assemble_jacobeans_1D(double* coordinates, dim_t numQuad,double* QuadWeights,
32                               dim_t numShape, dim_t numElements,index_t* nodes,
33                               double* DSDv, dim_t numTest,double* DTDv,
34                               double* dTdX, double* volume, index_t* element_id) {
35         #define DIM 1
36         #define LOCALDIM 1
37         register int e,q,s;
38         char error_msg[LenErrorMsg_MAX];
39         #pragma omp parallel
40         {
41           register double D,invD;
42           double X0[numShape];
43           #pragma omp for private(e,q,s) schedule(static)
44           for(e=0;e<numElements;e++){
45               for (s=0;s<numShape; s++) X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numShape)],DIM)];
46               for (q=0;q<numQuad;q++) {
47                  D=0;
48                  for (s=0;s<numShape; s++) D+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
49                  if (D==0.) {
50                      sprintf(error_msg,"Assemble_jacobeans_1D: element %d (id %d) has length zero.",e,element_id[e]);
51                      Finley_setError(ZERO_DIVISION_ERROR,error_msg);
52                  } else {
53                     invD=1./D;
54                     for (s=0;s<numTest; s++) dTdX[INDEX4(s,0,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*invD;
55             volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];
56                  }
57               }
58           }
59    
60  void Finley_Assemble_PDE(Finley_NodeFile* nodes,Finley_ElementFile* elements,Paso_SystemMatrix* S, escriptDataC* F,       }
61               escriptDataC* A, escriptDataC* B, escriptDataC* C, escriptDataC* D, escriptDataC* X, escriptDataC* Y ) {       #undef DIM
62         #undef LOCALDIM
63    
64    char error_msg[LenErrorMsg_MAX];  }
65    double *EM_S=NULL,*EM_F=NULL,*V=NULL,*dVdv=NULL,*dSdV=NULL,*Vol=NULL,*dvdV=NULL;  /**************************************************************/
66    double time0;  /*                                                            */
67    dim_t dimensions[ESCRIPT_MAX_DATA_RANK],e,q;  /*  Jacobean 2D                                               */
68    Assemble_Parameters p;  /*                                                            */
69    index_t *index_row=NULL,*index_col=NULL,color;  void Assemble_jacobeans_2D(double* coordinates, dim_t numQuad,double* QuadWeights,
70    Finley_resetError();                             dim_t numShape, dim_t numElements,index_t* nodes,
71                               double* DSDv, dim_t numTest,double* DTDv,
72    if (nodes==NULL || elements==NULL) return;                             double* dTdX, double* volume, index_t* element_id) {
73    if (S==NULL && isEmpty(F)) return;       #define DIM 2
74         #define LOCALDIM 2
75    /* set all parameters in p*/       register int e,q,s;
76    Assemble_getAssembleParameters(nodes,elements,S,F,&p);       char error_msg[LenErrorMsg_MAX];
77    if (! Finley_noError()) return;       #pragma omp parallel
78         {
79    /*  this function assumes NS=NN */         register double dXdv00,dXdv10,dXdv01,dXdv11,
80    if (p.NN!=p.NS) {                         dvdX00,dvdX10,dvdX01,dvdX11, D,invD;
81      Finley_setError(SYSTEM_ERROR,"Finley_Assemble_PDE: for Finley_Assemble_PDE numNodes and numShapes have to be identical.");         double X0[numShape], X1[numShape];
82      return;         #pragma omp for private(e,q,s) schedule(static)
83    }         for(e=0;e<numElements;e++){
84    if (p.numDim!=p.numElementDim) {             for (s=0;s<numShape; s++) {
85      Finley_setError(SYSTEM_ERROR,"Finley_Assemble_PDE: Finley_Assemble_PDE accepts volume elements only.");               X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numShape)],DIM)];
86      return;               X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numShape)],DIM)];
87    }             }
88    /*  get a functionspace */             for (q=0;q<numQuad;q++) {
89    type_t funcspace=UNKNOWN;                dXdv00=0;
90    updateFunctionSpaceType(funcspace,A);                dXdv10=0;
91    updateFunctionSpaceType(funcspace,B);                dXdv01=0;
92    updateFunctionSpaceType(funcspace,C);                dXdv11=0;
93    updateFunctionSpaceType(funcspace,D);                for (s=0;s<numShape; s++) {
94    updateFunctionSpaceType(funcspace,X);                   dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
95    updateFunctionSpaceType(funcspace,Y);                   dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
96    if (funcspace==UNKNOWN) return; /* all  data are empty */                   dXdv01+=X0[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
97                     dXdv11+=X1[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
98    /* check if all function spaces are the same */                }
99                  D  =  dXdv00*dXdv11 - dXdv01*dXdv10;
100    if (! functionSpaceTypeEqual(funcspace,A) ) {                if (D==0.) {
101          Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient A");                    sprintf(error_msg,"Assemble_jacobeans_2D: element %d (id %d) has area zero.",e,element_id[e]);
102    }                    Finley_setError(ZERO_DIVISION_ERROR,error_msg);
103    if (! functionSpaceTypeEqual(funcspace,B) ) {                } else {
104          Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient B");                   invD=1./D;
105    }                   dvdX00= dXdv11*invD;
106    if (! functionSpaceTypeEqual(funcspace,C) ) {                   dvdX10=-dXdv10*invD;
107          Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient C");                   dvdX01=-dXdv01*invD;
108    }                   dvdX11= dXdv00*invD;
109    if (! functionSpaceTypeEqual(funcspace,D) ) {  
110          Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient D");                   for (s=0;s<numTest; s++) {
111    }                     dTdX[INDEX4(s,0,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX10;
112    if (! functionSpaceTypeEqual(funcspace,X) ) {                     dTdX[INDEX4(s,1,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX11;
113          Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient X");                   }
114    }           volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];
115    if (! functionSpaceTypeEqual(funcspace,Y) ) {                }
116          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);  
   }  
   
   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);  
   }  
   
   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);  
   }  
   
   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);  
   }  
   
   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);  
   }  
   
   /*  check the dimensions: */  
     
   if (p.numEqu==1 && p.numComp==1) {  
     if (!isEmpty(A)) {  
       dimensions[0]=p.numDim;  
       dimensions[1]=p.numDim;  
       if (!isDataPointShapeEqual(A,2,dimensions)) {  
           sprintf(error_msg,"Finley_Assemble_PDE: coefficient A: illegal shape, expected shape (%d,%d)",dimensions[0],dimensions[1]);  
           Finley_setError(TYPE_ERROR,error_msg);  
       }  
     }  
     if (!isEmpty(B)) {  
        dimensions[0]=p.numDim;  
        if (!isDataPointShapeEqual(B,1,dimensions)) {  
           sprintf(error_msg,"Finley_Assemble_PDE: coefficient B: illegal shape (%d,)",dimensions[0]);  
           Finley_setError(TYPE_ERROR,error_msg);  
        }  
     }  
     if (!isEmpty(C)) {  
        dimensions[0]=p.numDim;  
        if (!isDataPointShapeEqual(C,1,dimensions)) {  
           sprintf(error_msg,"Finley_Assemble_PDE: coefficient C, expected shape (%d,)",dimensions[0]);  
           Finley_setError(TYPE_ERROR,error_msg);  
        }  
     }  
     if (!isEmpty(D)) {  
        if (!isDataPointShapeEqual(D,0,dimensions)) {  
           Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficient D, rank 0 expected.");  
117         }         }
118      }  
119      if (!isEmpty(X)) {       }
120         dimensions[0]=p.numDim;       #undef DIM
121         if (!isDataPointShapeEqual(X,1,dimensions)) {       #undef LOCALDIM
122            sprintf(error_msg,"Finley_Assemble_PDE: coefficient X, expected shape (%d,",dimensions[0]);  }
123            Finley_setError(TYPE_ERROR,error_msg);  /**************************************************************/
124    /*                                                            */
125    /*  Jacobean 3D                                               */
126    /*                                                            */
127    void Assemble_jacobeans_3D(double* coordinates, dim_t numQuad,double* QuadWeights,
128                               dim_t numShape, dim_t numElements,index_t* nodes,
129                               double* DSDv, dim_t numTest,double* DTDv,
130                               double* dTdX, double* volume, index_t* element_id) {
131         #define DIM 3
132         #define LOCALDIM 3
133         register int e,q,s;
134         char error_msg[LenErrorMsg_MAX];
135         #pragma omp parallel
136         {
137           register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,dXdv02,dXdv12,dXdv22,
138                           dvdX00,dvdX10,dvdX20,dvdX01,dvdX11,dvdX21,dvdX02,dvdX12,dvdX22, D,invD;
139           double X0[numShape], X1[numShape], X2[numShape];
140           #pragma omp for private(e,q,s) schedule(static)
141           for(e=0;e<numElements;e++){
142               for (s=0;s<numShape; s++) {
143                 X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numShape)],DIM)];
144                 X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numShape)],DIM)];
145                 X2[s]=coordinates[INDEX2(2,nodes[INDEX2(s,e,numShape)],DIM)];
146               }
147               for (q=0;q<numQuad;q++) {
148                  dXdv00=0;
149                  dXdv10=0;
150                  dXdv20=0;
151                  dXdv01=0;
152                  dXdv11=0;
153                  dXdv21=0;
154                  dXdv02=0;
155                  dXdv12=0;
156                  dXdv22=0;
157                  for (s=0;s<numShape; s++) {
158                     dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
159                     dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
160                     dXdv20+=X2[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
161                     dXdv01+=X0[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
162                     dXdv11+=X1[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
163                     dXdv21+=X2[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
164                     dXdv02+=X0[s]*DSDv[INDEX3(s,2,q,numShape,LOCALDIM)];
165                     dXdv12+=X1[s]*DSDv[INDEX3(s,2,q,numShape,LOCALDIM)];
166                     dXdv22+=X2[s]*DSDv[INDEX3(s,2,q,numShape,LOCALDIM)];
167                  }
168                  D  =  dXdv00*(dXdv11*dXdv22-dXdv12*dXdv21)+ dXdv01*(dXdv20*dXdv12-dXdv10*dXdv22)+dXdv02*(dXdv10*dXdv21-dXdv20*dXdv11);
169                  if (D==0.) {
170                      sprintf(error_msg,"Assemble_jacobeans_3D: element %d (id %d) has volume zero.",e,element_id[e]);
171                      Finley_setError(ZERO_DIVISION_ERROR,error_msg);
172                  } else {
173                     invD=1./D;
174                     dvdX00=(dXdv11*dXdv22-dXdv12*dXdv21)*invD;
175                     dvdX10=(dXdv20*dXdv12-dXdv10*dXdv22)*invD;
176                     dvdX20=(dXdv10*dXdv21-dXdv20*dXdv11)*invD;
177                     dvdX01=(dXdv02*dXdv21-dXdv01*dXdv22)*invD;
178                     dvdX11=(dXdv00*dXdv22-dXdv20*dXdv02)*invD;
179                     dvdX21=(dXdv01*dXdv20-dXdv00*dXdv21)*invD;
180                     dvdX02=(dXdv01*dXdv12-dXdv02*dXdv11)*invD;
181                     dvdX12=(dXdv02*dXdv10-dXdv00*dXdv12)*invD;
182                     dvdX22=(dXdv00*dXdv11-dXdv01*dXdv10)*invD;
183    
184                     for (s=0;s<numTest; s++) {
185                       dTdX[INDEX4(s,0,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX10+DTDv[INDEX3(s,2,q,numShape,LOCALDIM)]*dvdX20;
186                       dTdX[INDEX4(s,1,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX11+DTDv[INDEX3(s,2,q,numShape,LOCALDIM)]*dvdX21;
187                       dTdX[INDEX4(s,2,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX02+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX12+DTDv[INDEX3(s,2,q,numShape,LOCALDIM)]*dvdX22;
188                     }
189                 volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];
190                  }
191               }
192         }         }
193      }  
194      if (!isEmpty(Y)) {       }
195         if (!isDataPointShapeEqual(Y,0,dimensions)) {       #undef DIM
196            Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficient Y, rank 0 expected.");       #undef LOCALDIM
197    }
198    /**************************************************************/
199    /*                                                            */
200    /*  Jacobean 2D manifold in 3D                                */
201    /*                                                            */
202    void Assemble_jacobeans_3D_M2D(double* coordinates, dim_t numQuad,double* QuadWeights,
203                                   dim_t numShape, dim_t numElements,index_t* nodes,
204                                   double* DSDv, dim_t numTest,double* DTDv,
205                                   double* dTdX, double* volume, index_t* element_id) {
206         #define DIM 3
207         #define LOCALDIM 2
208         register int e,q,s;
209         char error_msg[LenErrorMsg_MAX];
210         #pragma omp parallel
211         {
212           register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,m00,m01,m11,
213                           dvdX00,dvdX01,dvdX02,dvdX10,dvdX11,dvdX12,D,invD;
214           double X0[numShape], X1[numShape], X2[numShape];
215           #pragma omp for private(e,q,s) schedule(static)
216           for(e=0;e<numElements;e++){
217               for (s=0;s<numShape; s++) {
218                 X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numShape)],DIM)];
219                 X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numShape)],DIM)];
220                 X2[s]=coordinates[INDEX2(2,nodes[INDEX2(s,e,numShape)],DIM)];
221               }
222               for (q=0;q<numQuad;q++) {
223                  dXdv00=0;
224                  dXdv10=0;
225                  dXdv20=0;
226                  dXdv01=0;
227                  dXdv11=0;
228                  dXdv21=0;
229                  for (s=0;s<numShape; s++) {
230                     dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
231                     dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
232                     dXdv20+=X2[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
233                     dXdv01+=X0[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
234                     dXdv11+=X1[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
235                     dXdv21+=X2[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
236                  }
237                  m00=dXdv00*dXdv00+dXdv10*dXdv10+dXdv20*dXdv20;
238                  m01=dXdv00*dXdv01+dXdv10*dXdv11+dXdv20*dXdv21;
239                  m11=dXdv01*dXdv01+dXdv11*dXdv11+dXdv21*dXdv21;
240                  D=m00*m11-m01*m01;
241                  if (D==0.) {
242                      sprintf(error_msg,"Assemble_jacobeans_3D_M2D: element %d (id %d) has area zero.",e,element_id[e]);
243                      Finley_setError(ZERO_DIVISION_ERROR,error_msg);
244                  } else {
245                     invD=1./D;
246                     dvdX00=( m00*dXdv00-m01*dXdv01)*invD;
247                     dvdX01=( m00*dXdv10-m01*dXdv11)*invD;
248                     dvdX02=( m00*dXdv20-m01*dXdv21)*invD;
249                     dvdX10=(-m01*dXdv00+m11*dXdv01)*invD;
250                     dvdX11=(-m01*dXdv10+m11*dXdv11)*invD;
251                     dvdX12=(-m01*dXdv20+m11*dXdv21)*invD;
252                     for (s=0;s<numTest; s++) {
253                       dTdX[INDEX4(s,0,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX10;
254                       dTdX[INDEX4(s,1,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX11;
255                       dTdX[INDEX4(s,2,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX02+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX12;
256                     }
257                 volume[INDEX2(q,e,numQuad)]=sqrt(D)*QuadWeights[q];
258                  }
259               }
260         }         }
261      }  
262    } else {       }
263      if (!isEmpty(A)) {       #undef DIM
264        dimensions[0]=p.numEqu;       #undef LOCALDIM
265        dimensions[1]=p.numDim;  }
266        dimensions[2]=p.numComp;  /**************************************************************/
267        dimensions[3]=p.numDim;  /*                                                            */
268        if (!isDataPointShapeEqual(A,4,dimensions)) {  /*  Jacobean 1D manifold in 3D                                */
269            sprintf(error_msg,"Finley_Assemble_PDE: coefficient A, expected shape (%d,%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2],dimensions[3]);  /*                                                            */
270            Finley_setError(TYPE_ERROR,error_msg);  void Assemble_jacobeans_3D_M1D(double* coordinates, dim_t numQuad,double* QuadWeights,
271        }                                 dim_t numShape, dim_t numElements,index_t* nodes,
272      }                                 double* DSDv, dim_t numTest,double* DTDv,
273      if (!isEmpty(B)) {                                 double* dTdX, double* volume, index_t* element_id) {
274        dimensions[0]=p.numEqu;       #define DIM 3
275        dimensions[1]=p.numDim;       #define LOCALDIM 1
276        dimensions[2]=p.numComp;       register int e,q,s;
277        if (!isDataPointShapeEqual(B,3,dimensions)) {       char error_msg[LenErrorMsg_MAX];
278            sprintf(error_msg,"Finley_Assemble_PDE: coefficient B, expected shape (%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2]);       #pragma omp parallel
           Finley_setError(TYPE_ERROR,error_msg);  
       }  
     }  
     if (!isEmpty(C)) {  
       dimensions[0]=p.numEqu;  
       dimensions[1]=p.numComp;  
       dimensions[2]=p.numDim;  
       if (!isDataPointShapeEqual(C,3,dimensions)) {  
           sprintf(error_msg,"Finley_Assemble_PDE: coefficient C, expected shape (%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2]);  
           Finley_setError(TYPE_ERROR,error_msg);  
       }  
     }  
     if (!isEmpty(D)) {  
       dimensions[0]=p.numEqu;  
       dimensions[1]=p.numComp;  
       if (!isDataPointShapeEqual(D,2,dimensions)) {  
           sprintf(error_msg,"Finley_Assemble_PDE: coefficient D, expected shape (%d,%d)",dimensions[0],dimensions[1]);  
           Finley_setError(TYPE_ERROR,error_msg);  
       }  
     }  
     if (!isEmpty(X)) {  
       dimensions[0]=p.numEqu;  
       dimensions[1]=p.numDim;  
       if (!isDataPointShapeEqual(X,2,dimensions)) {  
           sprintf(error_msg,"Finley_Assemble_PDE: coefficient X, expected shape (%d,%d)",dimensions[0],dimensions[1]);  
           Finley_setError(TYPE_ERROR,error_msg);  
       }  
     }  
     if (!isEmpty(Y)) {  
       dimensions[0]=p.numEqu;  
       if (!isDataPointShapeEqual(Y,1,dimensions)) {  
           sprintf(error_msg,"Finley_Assemble_PDE: coefficient Y, expected shape (%d,)",dimensions[0]);  
           Finley_setError(TYPE_ERROR,error_msg);  
       }  
     }  
   }  
   
   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)  
279       {       {
280           EM_S=EM_F=V=dVdv=dSdV=Vol=dvdV=NULL;         register double dXdv00,dXdv10,dXdv20,dvdX00,dvdX01,dvdX02,D,invD;
281           index_row=index_col=NULL;         double X0[numShape], X1[numShape], X2[numShape];
282           #pragma omp for private(e,q,s) schedule(static)
283           for(e=0;e<numElements;e++){
284               for (s=0;s<numShape; s++) {
285                 X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numShape)],DIM)];
286                 X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numShape)],DIM)];
287                 X2[s]=coordinates[INDEX2(2,nodes[INDEX2(s,e,numShape)],DIM)];
288               }
289               for (q=0;q<numQuad;q++) {
290                  dXdv00=0;
291                  dXdv10=0;
292                  dXdv20=0;
293                  for (s=0;s<numShape; s++) {
294                     dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
295                     dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
296                     dXdv20+=X2[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
297                  }
298                  D=dXdv00*dXdv00+dXdv10*dXdv10+dXdv20*dXdv20;
299                  if (D==0.) {
300                      sprintf(error_msg,"Assemble_jacobeans_3D_M1D: element %d (id %d) has length zero.",e,element_id[e]);
301                      Finley_setError(ZERO_DIVISION_ERROR,error_msg);
302                  } else {
303                     invD=1./D;
304                     dvdX00=dXdv00*invD;
305                     dvdX01=dXdv10*invD;
306                     dvdX02=dXdv20*invD;
307                     for (s=0;s<numTest; s++) {
308                       dTdX[INDEX4(s,0,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX00;
309                       dTdX[INDEX4(s,1,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX01;
310                       dTdX[INDEX4(s,2,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX02;
311                     }
312                 volume[INDEX2(q,e,numQuad)]=sqrt(D)*QuadWeights[q];
313                  }
314               }
315           }
316    
317           /* allocate work arrays: */       }
318         #undef DIM
319         #undef LOCALDIM
320    }
321    /**************************************************************/
322    /*                                                            */
323    /*  Jacobean 1D manifold in 2D                                */
324    /*                                                            */
325    void Assemble_jacobeans_2D_M1D(double* coordinates, dim_t numQuad,double* QuadWeights,
326                                   dim_t numShape, dim_t numElements,index_t* nodes,
327                                   double* DSDv, dim_t numTest,double* DTDv,
328                                   double* dTdX, double* volume, index_t* element_id) {
329         #define DIM 2
330         #define LOCALDIM 1
331         register int e,q,s;
332         char error_msg[LenErrorMsg_MAX];
333         #pragma omp parallel
334         {
335           register double dXdv00,dXdv10,dvdX00,dvdX01,D,invD;
336           double X0[numShape], X1[numShape];
337           #pragma omp for private(e,q,s) schedule(static)
338           for(e=0;e<numElements;e++){
339               for (s=0;s<numShape; s++) {
340                 X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numShape)],DIM)];
341                 X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numShape)],DIM)];
342               }
343               for (q=0;q<numQuad;q++) {
344                  dXdv00=0;
345                  dXdv10=0;
346                  for (s=0;s<numShape; s++) {
347                     dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
348                     dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
349                  }
350                  D=dXdv00*dXdv00+dXdv10*dXdv10;
351                  if (D==0.) {
352                      sprintf(error_msg,"Assemble_jacobeans_2D_M1D: element %d (id %d) has length zero.",e,element_id[e]);
353                      Finley_setError(ZERO_DIVISION_ERROR,error_msg);
354                  } else {
355                     invD=1./D;
356                     dvdX00=dXdv00*invD;
357                     dvdX01=dXdv10*invD;
358                     for (s=0;s<numTest; s++) {
359                       dTdX[INDEX4(s,0,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX00;
360                       dTdX[INDEX4(s,1,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX01;
361                     }
362                 volume[INDEX2(q,e,numQuad)]=sqrt(D)*QuadWeights[q];
363                  }
364               }
365           }
366    
          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: */  
 #ifndef PASO_MPI  
            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) {  
 #else  
            for(e=0;e<elements->numElements;e++){  
 #endif    
 //============================  
                   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(p.NN,&(elements->Nodes[INDEX2(0,e,p.NN)]),p.numDim,nodes->Coordinates,V);  
                   /*  calculate dVdv(i,j,q)=V(i,k)*DSDv(k,j,q) */  
           Finley_Util_SmallMatMult(p.numDim,p.numDim*p.numQuad,dVdv,p.NS,V,p.referenceElement->dSdv);  
                   /*  dvdV=invert(dVdv) inplace */  
           Finley_Util_InvertSmallMat(p.numQuad,p.numDim,dVdv,dvdV,Vol);  
                   /*  calculate dSdV=DSDv*DvDV */  
           Finley_Util_SmallMatSetMult(p.numQuad,p.NS_row,p.numDim,dSdV,p.numDim,p.referenceElement_row->dSdv,dvdV);  
                   /*  scale volume: */  
           for (q=0;q<p.numQuad;q++) Vol[q]=ABS(Vol[q]*p.referenceElement->QuadWeights[q]);  
 //============================  
       
                    /*   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  */  
 #ifndef PASO_MPI  
                        Finley_Util_AddScatter(p.NN_row,index_row,p.numEqu,EM_F,getSampleData(F,0));  
 #else  
                        Finley_Util_AddScatter_upperBound(p.NN_row,index_row,p.numEqu,EM_F,getSampleData(F,0),p.degreeOfFreedomUpperBound);  
 #endif                                        
                     }  
                 }  
               }  
 #ifndef PASO_MPI  
             }  
          }  
 #endif  
          /* 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);  
367       }       }
368       #ifdef Finley_TRACE       #undef DIM
369       printf("timing: assemblage PDE: %.4e sec\n",Finley_timer()-time0);       #undef LOCALDIM
      #endif  
   }  
370  }  }
371  /*  /*
372   * $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  
  *  
  *  
373   *   *
374   */   */

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

  ViewVC Help
Powered by ViewVC 1.1.26