Charm++: Charm Source Code Documentation

00001 
00008 #include <charm.h> /* for CkAbort */
00009 #include <vector>
00010 #include <list>
00011 using std::list;
00012 using std::vector;
00013 
00014 #include <GenericElement.h>
00015 
00016 /* Numerical fudge factors */
00017 #ifndef TOL
00018 #define TOL 1e-9
00019 #endif
00021 const double LTOL = 0.0 - TOL;
00022 const double HTOL = 1.0 + TOL;
00023 
00024 /************* Element Shape Functions **********/
00025 
00026 // Populate SF with the shape function of the element
00027 void 
00028 GenericElement::shape_func(const CVector &nc,double SF[]) const
00029 {
00030   switch (_size) {
00031     // Tets are 4 nodes or 10
00032   case 4: {
00033     SF[0] = 1. - nc.x - nc.y - nc.z;
00034     SF[1] = nc.x;
00035     SF[2] = nc.y;
00036     SF[3] = nc.z;
00037     break;
00038   }
00039   case 10: {
00040     const double xi = nc.x;
00041     const double eta = nc.y;
00042     const double zeta = nc.z;
00043     const double alpha = (1. - xi - eta - zeta);
00044     SF[0] = alpha*(1. - 2.*(xi + eta + zeta));
00045     SF[1] = xi *(2.* xi - 1.);
00046     SF[2] = eta *(2. * eta - 1.);
00047     SF[3] = zeta *(2. * zeta - 1.);
00048     SF[4] = 4.* xi * alpha;
00049     SF[5] = 4.* eta * alpha;
00050     SF[6] = 4.* zeta * alpha;
00051     SF[7] = 4. * xi * eta;
00052     SF[8] = 4. * eta * zeta;
00053     SF[9] = 4. * xi * zeta;
00054     break;
00055   }
00056     // Hex's are 8 nodes or 20
00057   case 8: {
00058     const double xi = nc.x;
00059     const double xi_minus = 1. - xi;
00060     const double eta = nc.y;
00061     const double eta_minus = 1. - eta;
00062     const double zeta = nc.z;
00063     const double zeta_minus = 1. - zeta;
00064     SF[0] = xi_minus * eta_minus * zeta_minus;
00065     SF[1] = xi * eta_minus * zeta_minus;
00066     SF[2] = xi * eta * zeta_minus;
00067     SF[3] = xi_minus * eta * zeta_minus;
00068     SF[4] = xi_minus * eta_minus * zeta;
00069     SF[5] = xi * eta_minus * zeta;
00070     SF[6] = xi * eta * zeta;
00071     SF[7] = xi_minus * eta * zeta;
00072     break;
00073   }
00074   default: 
00075     CkAbort("GenericElement::shape_func:Error: unkown element type.");
00076   }
00077 }
00078 
00079 // Populate dSF with the derivative of the shape function for the element
00080 void
00081 GenericElement::dshape_func(const CVector &nc,double dSF[][3]) const
00082 {
00083   switch (_size) {
00084   case 4: {
00085     dSF[0][0] = -1;  dSF[0][1] = -1;  dSF[0][2] = -1;
00086     dSF[1][0] =  1;  dSF[1][1] =  0;  dSF[1][2] =  0;
00087     dSF[2][0] =  0;  dSF[2][1] =  1;  dSF[2][2] =  0;
00088     dSF[3][0] =  0;  dSF[3][1] =  0;  dSF[3][2] =  1;
00089     break;
00090   }
00091   case 10:{
00092     const double xi = nc.x;
00093     const double eta = nc.y;
00094     const double zeta = nc.z;
00095     const double alpha = (1. - xi - eta - zeta);
00096     dSF[0][0] = (4.*(xi+eta+zeta)-3.);      dSF[0][1] = dSF[0][0];                 dSF[0][2] = dSF[0][0];
00097     dSF[1][0] = 4.*xi - 1.;                 dSF[1][1] = 0;                         dSF[1][2] = 0;
00098     dSF[2][0] = 0;                          dSF[2][1] = 4.*eta - 1.;               dSF[2][2] = 0;
00099     dSF[3][0] = 0;                          dSF[3][1] = 0;                         dSF[3][2] = 4.*zeta - 1.;
00100     dSF[4][0] = 4.*(alpha - xi);            dSF[4][1] = -4.*xi;                    dSF[4][2] = -4.*xi;
00101     dSF[5][0] = -4.*eta;                    dSF[5][1] = 4.*(alpha - eta);          dSF[5][2] = -4.*eta;
00102     dSF[6][0] = -4.*zeta;                   dSF[6][1] = -4.*zeta;                  dSF[6][2] = 4.*(alpha - zeta);
00103     dSF[7][0] = 4.*eta;                     dSF[7][1] = 4.*xi;                     dSF[7][2] = 0;
00104     dSF[8][0] = 0;                          dSF[8][1] = 4.*zeta;                   dSF[8][2] = 4.*eta;
00105     dSF[9][0] = 4.*zeta;                    dSF[9][1] = 0;                         dSF[9][2] = 4.*xi;
00106     break;
00107   }    
00108   case 8: {
00109     const double xi = nc.x;
00110     const double xi_minus = 1. - xi;
00111     const double eta = nc.y;
00112     const double eta_minus = 1. - eta;
00113     const double zeta = nc.z;
00114     const double zeta_minus = 1. - zeta;
00115     dSF[0][0] = -1.*eta_minus*zeta_minus;  dSF[0][1] = -1.*xi_minus*zeta_minus;  dSF[0][2] = -1.*xi_minus*eta_minus;
00116     dSF[1][0] = eta_minus*zeta_minus;      dSF[1][1] = -1.*xi*zeta_minus;        dSF[1][2] = -1.*xi*eta_minus;
00117     dSF[2][0] = eta*zeta_minus;            dSF[2][1] = xi*zeta_minus;            dSF[2][2] = -1.*xi*eta;
00118     dSF[3][0] = -1.*eta*zeta_minus;        dSF[3][1] = xi_minus*zeta_minus;      dSF[3][2] = -1.*xi_minus*eta;
00119     dSF[4][0] = -1.*eta_minus*zeta;        dSF[4][1] = -1.*xi_minus*zeta;        dSF[4][2] = xi_minus*eta_minus;
00120     dSF[5][0] = eta_minus*zeta;            dSF[5][1] = -1.*xi*zeta;              dSF[5][2] = xi*eta_minus;
00121     dSF[6][0] = eta*zeta;                  dSF[6][1] = xi*zeta;                  dSF[6][2] = xi*eta;
00122     dSF[7][0] = -1.*eta*zeta;              dSF[7][1] = xi_minus * zeta;          dSF[7][2] = xi_minus*eta;
00123     break;
00124   }
00125   default:
00126     CkAbort("GenericElement::dshape_func:error: Unknown element type.");
00127   }
00128 }
00129 
00130 // Populate J with the jacobian for the point p having natural coordinates nc.  
00131 void
00132 GenericElement::jacobian(const CPoint p[],const CVector &nc,CVector J[]) const
00133 {
00134   switch(_size){
00135   case 4: {
00136     J[0] = p[1] - p[0];
00137     J[1] = p[2] - p[0];
00138     J[2] = p[3] - p[0];
00139     break;
00140   }
00141   case 10: {
00142     const double xi = nc.x;
00143     const double eta = nc.y;
00144     const double zeta = nc.z;
00145     const double alpha = (1. - xi - eta - zeta);
00146     CPoint P(p[0]*(4.*(xi+eta+zeta)-3.));
00147     J[0] = ((p[9]-p[6])*4.*zeta)+((p[7]-p[5])*4.*eta)+
00148       (p[4]*(4.*(alpha-xi))+p[1]*(4.*xi-1.)+P);
00149     J[1] = ((p[8]-p[6])*4.*zeta)+((p[7]-p[4])*4.*xi)+
00150       (p[5]*(4.*(alpha-eta))+p[2]*(4.*eta-1.)+P);
00151     J[2] = ((p[9]-p[4])*4.*xi)+((p[8]-p[5])*4.*eta)+
00152       (p[6]*(4.*(alpha-zeta))+p[3]*(4.*zeta-1.)+P);
00153     break;
00154   }
00155   case 8: {
00156     const double xi = nc.x;
00157     const double xi_minus = 1. - xi;
00158     const double eta = nc.y;
00159     const double eta_minus = 1. - eta;
00160     const double zeta = nc.z;
00161     const double zeta_minus = 1. - zeta;
00162     J[0] = ((p[6]-p[7])*eta*zeta)+((p[5]-p[4])*eta_minus*zeta)+
00163       ((p[2]-p[3])*eta*zeta_minus)+((p[1]-p[0])*eta_minus*zeta_minus);
00164     J[1] = ((p[7]-p[4])*xi_minus*zeta)+((p[6]-p[5])*xi*zeta)+
00165       ((p[3]-p[0])*xi_minus*zeta_minus)+((p[2]-p[1])*xi*zeta_minus);
00166     J[2] = ((p[7]-p[3])*xi_minus*eta)+((p[6]-p[2])*xi*eta)+
00167       ((p[5]-p[1])*xi*eta_minus)+((p[4]-p[0])*xi_minus*eta_minus);
00168     break;
00169   }
00170   default:
00171     CkAbort("GenericElement::jacobian:Error: Cannot handle this element size (yet).");
00172   }
00173 }
00174 
00176 void 
00177 GenericElement::interpolate_natural(int nValuesPerNode,
00178            const ConcreteElementNodeData &src, // Source element
00179            const CVector &nc,
00180            double *dest) const
00181 {
00182   const double xi = nc.x;
00183   const double eta = nc.y;
00184   const double zeta = nc.z;
00185   for (int i=0;i<nValuesPerNode;i++) // loop over data values
00186   {
00187     double f[maxSize]; // data for each of our nodes
00188     for (int n=0;n<_size;n++) f[n]=src.getNodeData(n)[i];
00189     switch(_size) {
00190     case 4: {
00191       dest[i] = f[0]+(((f[1]-f[0])*xi) + ((f[2]-f[0])*eta) + ((f[3] - f[0])*zeta));
00192       break;
00193     }
00194     case 10: {
00195       const double alpha = (1.-xi-eta-zeta);
00196       dest[i] = (alpha*(1.-2.*(xi+eta+zeta))*f[0] +
00197        xi*(2.*xi-1.)*f[1] +
00198        eta*(2.*eta-1.)*f[2] +
00199        zeta*(2.*zeta-1.)*f[3] +
00200        4.*xi*alpha*f[4] +
00201        4.*eta*alpha*f[5] +
00202        4.*zeta*alpha*f[6] +
00203        4.*xi*eta*f[7] +
00204        4.*eta*zeta*f[8] +
00205        4.*zeta*xi*f[9]);
00206       break;
00207     }
00208     case 8: {
00209       const double xi = nc.x;
00210       const double xi_minus = 1. - xi;
00211       const double eta = nc.y;
00212       const double eta_minus = 1. - eta;
00213       const double zeta = nc.z;
00214       const double zeta_minus = 1. - zeta;
00215       dest[i] = (xi_minus*eta_minus*zeta_minus*f[0] +
00216        xi*eta_minus*zeta_minus*f[1] +
00217        xi*eta*zeta_minus*f[2] +
00218        xi_minus*eta*zeta_minus*f[3] +
00219        xi_minus*eta_minus*zeta*f[4] +
00220        xi*eta_minus*zeta*f[5] +
00221        xi*eta*zeta*f[6] +
00222        xi_minus*eta*zeta*f[7]);
00223       break;
00224     }
00225     default:
00226       CkAbort("interpolate::error Cannot handle this element type (yet).");
00227     }
00228   }
00229 }
00230 
00232 void 
00233 Transpose(CVector matrix[])
00234 {
00235   CVector tpose[3];
00236   
00237   tpose[0]=CVector(matrix[0].x,matrix[1].x,matrix[2].x);
00238   tpose[1]=CVector(matrix[0].y,matrix[1].y,matrix[2].y);
00239   tpose[2]=CVector(matrix[0].z,matrix[1].z,matrix[2].z);
00240   matrix[0]=tpose[0];
00241   matrix[1]=tpose[1];
00242   matrix[2]=tpose[2];
00243 }
00244 
00245 // Evaluates the shape function and jacobian at a given point
00246 void
00247 GenericElement::shapef_jacobian_at(const CPoint &p,CVector &natc,
00248                    const ConcreteElement &e, // Source element
00249                    CVector &fvec,CVector fjac[]) const
00250 {  
00251   CPoint P[maxSize];
00252   double SF[maxSize];
00253   this->shape_func(natc,SF);
00254   fvec=-p;
00255   for(int i = 0;i < _size;i++){
00256     P[i]=e.getNodeLocation(i);
00257     fvec += SF[i]*P[i];
00258   }
00259   this->jacobian(P,natc,fjac);
00260   Transpose(fjac);
00261 }
00262 
00263 
00264 
00265 /******************* Point Location ********************/
00266 
00267 
00268 bool
00269 NewtonRaphson(CVector &natc,
00270           const GenericElement &el,
00271           const ConcreteElement &e, // Source element
00272           const CPoint &point);
00273 bool 
00274 LUDcmp(CVector a[], 
00275        int indx[]);
00276 
00277 void 
00278 LUBksb(CVector a[],
00279        int indx[], 
00280        CVector &b);
00281 
00282 
00283 // Return true and the natural coordinates if this point lies in this element.
00284 bool
00285 GenericElement::element_contains_point(const CPoint &p, //    Target Mesh point
00286                    const ConcreteElement &e, // Source element
00287                    CVector &natc) const // Returns Targ nat
00288 {
00289   // guess at natural coordinates of point: center of element
00290   if(_size == 4 || _size == 10)
00291     natc=CVector(.25,.25,.25);
00292   else if (_size == 8 || _size == 20)
00293     natc=CVector(.5,.5,.5);
00294   else{
00295     CkAbort("GenericElement::element_contains_point: Error: Cannot handle this element type. (yet)");
00296   }
00297   
00298   // Solve for the natural coordinates using non-linear newton-raphson 
00299   if(!NewtonRaphson(natc,*this,e,p)){
00300     CkAbort("GenericElement::global_find_point_in_mesh: error NewtonRaphson failed.");
00301   }
00302   
00303   // Make sure natural coordinates lie in unit cube:
00304   if(natc[0] >= LTOL && natc[0] <= HTOL &&
00305      natc[1] >= LTOL && natc[1] <= HTOL &&
00306      natc[2] >= LTOL && natc[2] <= HTOL)
00307   {
00308     if(_size == 4 || _size == 10){
00309       // Natural coordinates must sum to less than one.
00310       if((natc[0]+natc[1]+natc[2]) <= HTOL){
00311         return (true);
00312       }
00313     }
00314     else if(_size == 8 || _size == 20) { 
00315       // Everything in unit cube is OK.
00316       return(true);
00317     } 
00318   }
00319   return false; // Point's natural coordinates do not lie in unit cube
00320 }
00321 
00322 // Newton-Raphson method, customized for using the GeoPrimitives
00323 // and computational mesh constructs
00324 //
00325 // Parameters: 
00326 // double natc[] = initial guess at the natural coordinates 
00327 // unsigned int elnum = index of the mesh element to use
00328 // const ElementConnectivity &ec = The element connectivity for the mesh
00329 // const NodalCoordinates &nc = The nodal coordinates for the mesh
00330 // const CPoint &point = the point at which we wish a solution
00331 //
00332 bool
00333 NewtonRaphson(CVector &natc,
00334           const GenericElement &el,
00335           const ConcreteElement &e, // Source element
00336           const CPoint &point) 
00337 {
00338   int k,i;
00339   int ntrial = 4; // Number of iterations before giving up
00340   double errx,errf,d;
00341   int indx[3];
00342   CVector p;
00343   CVector fvec;
00344   CVector fjac[3];
00345   for (k=0;k<ntrial;k++) {
00346     el.shapef_jacobian_at(point,natc,e,fvec,fjac);
00347     errf=0.0;
00348     for (i=0;i<3;i++) 
00349       errf += fabs(fvec[i]);
00350     if (errf <= TOL)
00351       return (true);
00352     p = -1.0 * fvec;
00353     if(!LUDcmp(fjac,indx)){
00354       // cerr << "NewtonRaphson::error: LUDcmp failed." << endl;
00355       return(false);
00356     }
00357     LUBksb(fjac,indx,p);
00358     errx=0.0;
00359     for (i=0;i<3;i++) 
00360       errx += fabs(p[i]);
00361     natc += p;
00362     if (errx <= TOL)
00363       return (true);
00364   }
00365   // cerr << "NewtonRaphson::warning: reached maximum iterations" << endl;
00366   return (true);
00367 }
00368 
00369 
00370 // LU Decomp
00371 #define TINY 1.0e-20
00372 bool
00373 LUDcmp(CVector a[], int indx[])
00374 {
00375   int i,imax,j,k;
00376   double big,dum,sum,temp,d;
00377   CVector vv;
00378   
00379   
00380   
00381   for (i=0;i<3;i++) {
00382     big=0.0;
00383     for (j=0;j<3;j++)
00384       if ((temp=fabs(a[i][j])) > big) big=temp;
00385     if (big == 0.0){
00386       // cerr << "LUDcmp::error: Singular matrix" << endl;
00387       return(false);
00388     }
00389     vv[i]=1.0/big;
00390   }
00391   for (j=0;j<3;j++) {
00392     for (i=0;i<j;i++) {
00393       sum=a[i][j];
00394       for (k=0;k<i;k++) 
00395     sum -= a[i][k]*a[k][j];
00396       a[i][j]=sum;
00397     }
00398     big=0.0;
00399     for (i=j;i<3;i++) {
00400       sum=a[i][j];
00401       for (k=0;k<j;k++)
00402     sum -= a[i][k]*a[k][j];
00403       a[i][j]=sum;
00404       if ( (dum=vv[i]*fabs(sum)) >= big) {
00405     big=dum;
00406     imax=i;
00407       }
00408     }
00409     if (j != imax) {
00410       for (k=0;k<3;k++) {
00411     dum=a[imax][k];
00412     a[imax][k]=a[j][k];
00413     a[j][k]=dum;
00414       }
00415       vv[imax]=vv[j];
00416     }
00417     indx[j]=imax;
00418     if (a[j][j] == 0.0) a[j][j]=TINY;
00419     if (j != 3) {
00420       dum=1.0/(a[j][j]);
00421       for (i=j+1;i<3;i++) 
00422     a[i][j] *= dum;
00423     }
00424   }
00425   return (true);
00426 }
00427 #undef TINY
00428 
00429 // LU Back substitution
00430 void
00431 LUBksb(CVector a[],int indx[], CVector &b)
00432 {
00433   int i,ii=0,ip,j;
00434   double sum;
00435   
00436   for (i=0;i<3;i++) {
00437     ip=indx[i];
00438     sum=b[ip];
00439     b[ip]=b[i];
00440     if (ii){
00441       for (j=ii;j<=i-1;j++) 
00442     sum -= a[i][j]*b[j];
00443     }
00444     else if (sum) 
00445       ii=i;
00446     b[i]=sum;
00447   }
00448   for (i=2;i>=0;i--) {
00449     sum=b[i];
00450     for (j=i+1;j<3;j++) 
00451       sum -= a[i][j]*b[j];
00452     b[i]=sum/a[i][i];
00453   }
00454 }
00455
libs/ck-libs/datatransfer/GenericElement.C
