Changeset 383


Ignore:
Timestamp:
05/13/09 10:11:32 (16 years ago)
Author:
seroussi
Message:

we have now a 2d an 3d determinant for the tria

Location:
issm/trunk/src/c/objects
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/objects/Icefront.cpp

    r246 r383  
    899899                                complete_list[j][2]=0;
    900900                        }
    901                         tria->GetJacobianDeterminant(&J[i],&complete_list[0][0],l1l2l3_tria);
     901                        tria->GetJacobianDeterminant3d(&J[i],&complete_list[0][0],l1l2l3_tria);
    902902                }
    903903

  • issm/trunk/src/c/objects/Tria.cpp

    r358 r383  
    410410               
    411411                /* Get Jacobian determinant: */
    412                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     412                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    413413
    414414                /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant
    … …  
    640640
    641641                /* Get Jacobian determinant: */
    642                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     642                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    643643
    644644                /*Get L matrix: */
    … …  
    732732
    733733                /* Get Jacobian determinant: */
    734                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     734                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    735735               
    736736                DL_scalar=gauss_weight*Jdet;
    … …  
    882882
    883883                /* Get Jacobian determinant: */
    884                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     884                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    885885               
    886886                 /*Get nodal functions: */
    … …  
    10141014               
    10151015                /* Get Jacobian determinant: */
    1016                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     1016                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    10171017               
    10181018                 /*Get nodal functions: */
    … …  
    11941194
    11951195#undef __FUNCT__
    1196 #define __FUNCT__ "Tria::GetJacobianDeterminant"
    1197 void Tria::GetJacobianDeterminant(double*  Jdet, double* xyz_list,double* gauss_l1l2l3){
     1196#define __FUNCT__ "Tria::GetJacobianDeterminant2d"
     1197void Tria::GetJacobianDeterminant2d(double*  Jdet, double* xyz_list,double* gauss_l1l2l3){
     1198
     1199        /*The Jacobian determinant is constant over the element, discard the gaussian points.
     1200         * J is assumed to have been allocated of size NDOF2xNDOF2.*/
     1201
     1202        double x1,x2,x3,y1,y2,y3;
     1203       
     1204        x1=*(xyz_list+3*0+0);
     1205        y1=*(xyz_list+3*0+1);
     1206        x2=*(xyz_list+3*1+0);
     1207        y2=*(xyz_list+3*1+1);
     1208        x3=*(xyz_list+3*2+0);
     1209        y3=*(xyz_list+3*2+1);
     1210
     1211
     1212        *Jdet=sqrt(3.0)/6.0*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1));
     1213
     1214
     1215        if(Jdet<0){
     1216                printf("%s%s\n",__FUNCT__," error message: negative jacobian determinant!");
     1217        }
     1218       
     1219}
     1220
     1221#undef __FUNCT__
     1222#define __FUNCT__ "Tria::GetJacobianDeterminant3d"
     1223void Tria::GetJacobianDeterminant3d(double*  Jdet, double* xyz_list,double* gauss_l1l2l3){
    11981224
    11991225        /*The Jacobian determinant is constant over the element, discard the gaussian points.
    … …  
    16641690
    16651691                /* Get Jacobian determinant: */
    1666                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     1692                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    16671693                #ifdef _DEBUG_
    16681694                printf("Element id %i Jacobian determinant: %lf\n",TriaElementGetID(this),Jdet);
    … …  
    19121938
    19131939                /* Get Jacobian determinant: */
    1914                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     1940                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    19151941                #ifdef _DEBUG_
    19161942                printf("Element id %i Jacobian determinant: %lf\n",TriaElementGetID(this),Jdet);
    … …  
    20532079
    20542080                /* Get Jacobian determinant: */
    2055                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     2081                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    20562082               
    20572083                /* Get nodal functions value at gaussian point:*/
    … …  
    21792205
    21802206                /* Get Jacobian determinant: */
    2181                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     2207                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    21822208                #ifdef _DEBUG_
    21832209                printf("Element id %i Jacobian determinant: %lf\n",TriaElementGetID(this),Jdet);
    … …  
    22452271
    22462272
    2247 
    2248 #undef __FUNCT__
    2249 #define __FUNCT__ "Tria::CreateKMatrixDiagnosticBaseVert"
    2250 void  Tria::CreateKMatrixDiagnosticBaseVert(Mat Kgg,void* vinputs,int analysis_type){
    2251 
    2252         int i,j;
    2253 
    2254         /* node data: */
    2255         const int    numgrids=3;
    2256         const int    NDOF1=1;
    2257         const int    numdofs=NDOF1*numgrids;
    2258         double       xyz_list[numgrids][3];
    2259         int          doflist[numdofs];
    2260         int          numberofdofspernode;
    2261        
    2262         /* gaussian points: */
    2263         int     num_gauss,ig;
    2264         double* first_gauss_area_coord  =  NULL;
    2265         double* second_gauss_area_coord =  NULL;
    2266         double* third_gauss_area_coord  =  NULL;
    2267         double* gauss_weights           =  NULL;
    2268         double  gauss_weight;
    2269         double  gauss_l1l2l3[3];
    2270 
    2271 
    2272         /* matrices: */
    2273         double L[1][3];
    2274         double DL_scalar;
    2275 
    2276         double Ke_gg[3][3]; //stiffness matrix
    2277         double Ke_gg_gaussian[3][3]; //stiffness matrix evaluated at the gaussian point.
    2278         double Jdet;
    2279 
    2280         ParameterInputs* inputs=NULL;
    2281 
    2282         /*recover pointers: */
    2283         inputs=(ParameterInputs*)vinputs;
    2284 
    2285         /* Get node coordinates and dof list: */
    2286         GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
    2287         GetDofList(&doflist[0],&numberofdofspernode);
    2288 
    2289         /* Set Ke_gg to 0: */
    2290         for(i=0;i<numdofs;i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]=0.0;
    2291 
    2292         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    2293         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    2294 
    2295         /* Start  looping on the number of gaussian points: */
    2296         for (ig=0; ig<num_gauss; ig++){
    2297                 /*Pick up the gaussian point: */
    2298                 gauss_weight=*(gauss_weights+ig);
    2299                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    2300                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    2301                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
    2302 
    2303                 /* Get Jacobian determinant: */
    2304                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    2305        
    2306                 //Get L matrix if viscous basal drag present:
    2307                 GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,NDOF1);
    2308 
    2309                 DL_scalar=gauss_weight*Jdet;
    2310                
    2311                 /*  Do the triple producte tL*D*L: */
    2312                 TripleMultiply( &L[0][0],1,3,1,
    2313                                         &DL_scalar,1,1,0,
    2314                                         &L[0][0],1,3,0,
    2315                                         &Ke_gg_gaussian[0][0],0);
    2316 
    2317                
    2318                 /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
    2319                 for( i=0; i<numdofs; i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
    2320         } //for (ig=0; ig<num_gauss; ig++
    2321 
    2322         /*Add Ke_gg to global matrix Kgg: */
    2323         MatSetValues(Kgg,numdofs,doflist,numdofs,doflist,(const double*)Ke_gg,ADD_VALUES);
    2324                
    2325         cleanup_and_return:
    2326         xfree((void**)&first_gauss_area_coord);
    2327         xfree((void**)&second_gauss_area_coord);
    2328         xfree((void**)&third_gauss_area_coord);
    2329         xfree((void**)&gauss_weights);
    2330 }
    2331 
    23322273#undef __FUNCT__
    23332274#define __FUNCT__ "Tria::CreateKMatrixDiagnosticSurfaceVert"
    … …  
    24262367
    24272368                /* Get Jacobian determinant: */
    2428                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     2369                GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    24292370       
    24302371                //Get L matrix if viscous basal drag present:
    … …  
    25472488
    25482489                /* Get Jacobian determinant: */
    2549                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     2490                GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    25502491               
    25512492                //Get L matrix if viscous basal drag present:

  • issm/trunk/src/c/objects/Tria.h

    r358 r383  
    7272                void  CreateKMatrixDiagnosticHoriz(Mat Kgg,void* inputs,int analysis_type);
    7373                void  CreateKMatrixDiagnosticHorizFriction(Mat Kgg,void* inputs,int analysis_type);
    74                 void  CreateKMatrixDiagnosticBaseVert(Mat Kgg,void* inputs,int analysis_type);
    7574                void  CreateKMatrixDiagnosticSurfaceVert(Mat Kgg,void* inputs,int analysis_type);
    7675                void  CreateKMatrixSlopeCompute(Mat Kgg,void* vinputs,int analysis_type);
    … …  
    8685                void  GetNodalFunctionsDerivativesParams(double* dl1dl2dl3,double* gauss_l1l2l3);
    8786                void  GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss_l1l2l3);
    88                 void  GetJacobian(double* J, double* xyz_list,double* gauss_l1l2l3);
     87                void  GetJacobian2d(double* J, double* xyz_list,double* gauss_l1l2l3);
     88                void  GetJacobian3d(double* J, double* xyz_list,double* gauss_l1l2l3);
    8989                void  Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type);
    9090                void  Gradj(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,char* control_type);
Note: See TracChangeset for help on using the changeset viewer.