Changeset 5625


Ignore:
Timestamp:
08/30/10 17:01:51 (15 years ago)
Author:
Mathieu Morlighem
Message:

Added GaussTria for slope only for now

Location:
issm/trunk/src/c
Files:
3 added
5 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/Makefile.am

    r5595 r5625  
    7777                                        ./objects/Bamg/Mesh.cpp\
    7878                                        ./objects/Bamg/Mesh.h\
     79                                        ./objects/Gauss/GaussTria.h\
     80                                        ./objects/Gauss/GaussTria.cpp\
    7981                                        ./objects/Update.h\
    8082                                        ./objects/Element.h\
    … …  
    630632                                        ./objects/Bamg/Mesh.h\
    631633                                        ./objects/Bamg/Mesh.cpp\
     634                                        ./objects/Gauss/GaussTria.h\
     635                                        ./objects/Gauss/GaussTria.cpp\
    632636                                        ./objects/Update.h\
    633637                                        ./objects/Element.h\

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

    r5624 r5625  
    45804580void  Tria::CreateKMatrixSlope(Mat Kgg){
    45814581
    4582         /* local declarations */
    4583         int             i,j;
    4584 
    4585         /* node data: */
    4586         const int    numgrids=3;
     4582        /*constants: */
     4583        const int    numnodes=3;
    45874584        const int    NDOF1=1;
    4588         const int    numdof=NDOF1*numgrids;
    4589         double       xyz_list[numgrids][3];
    4590         int*         doflist=NULL;
    4591        
    4592         /* gaussian points: */
    4593         int     num_gauss,ig;
    4594         double* first_gauss_area_coord  =  NULL;
    4595         double* second_gauss_area_coord =  NULL;
    4596         double* third_gauss_area_coord  =  NULL;
    4597         double* gauss_weights           =  NULL;
    4598         double  gauss_weight;
    4599         double  gauss_l1l2l3[3];
    4600 
    4601         /* matrices: */
    4602         double L[1][3];
    4603         double DL_scalar;
    4604 
    4605         /* local element matrices: */
    4606         double Ke_gg[numdof][numdof]={0.0}; //local element stiffness matrix
    4607         double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
    4608        
    4609         double Jdet;
    4610 
    4611         /* Get node coordinates and dof list: */
    4612         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     4585        const int    numdof=NDOF1*numnodes;
     4586
     4587        /* Intermediaries */
     4588        int        i,j,ig;
     4589        double     DL_scalar,Jdet;
     4590        double     xyz_list[numnodes][3];
     4591        double     L[1][3];
     4592        double     Ke[numdof][numdof] = {0.0};
     4593        double     Ke_g[numdof][numdof];
     4594        GaussTria *gauss = NULL;
     4595        int       *doflist = NULL;
     4596
     4597        GetVerticesCoordinates(&xyz_list[0][0], nodes, numnodes);
    46134598        GetDofList(&doflist);
    46144599
    4615         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    4616         GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    4617 
    4618         /* Start  looping on the number of gaussian points: */
    4619         for (ig=0; ig<num_gauss; ig++){
    4620                 /*Pick up the gaussian point: */
    4621                 gauss_weight=*(gauss_weights+ig);
    4622                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    4623                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    4624                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
    4625 
     4600        /* Start looping on the number of gaussian points: */
     4601        gauss=new GaussTria(2);
     4602        for (ig=gauss->begin();ig<gauss->end();ig++){
     4603
     4604                gauss->GaussPoint(ig);
    46264605               
    4627                 /*Get L matrix: */
    4628                 GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,NDOF1);
    4629 
    4630                 /* Get Jacobian determinant: */
    4631                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    4632                
    4633                 DL_scalar=gauss_weight*Jdet;
    4634 
    4635                 /*  Do the triple producte tL*D*L: */
    4636                 TripleMultiply( &L[0][0],1,3,1,
     4606                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
     4607                DL_scalar=gauss->weight*Jdet;
     4608
     4609                GetL(&L[0][0], &xyz_list[0][0], gauss,NDOF1);
     4610
     4611                TripleMultiply(&L[0][0],1,3,1,
    46374612                                        &DL_scalar,1,1,0,
    46384613                                        &L[0][0],1,3,0,
    4639                                         &Ke_gg_gaussian[0][0],0);
    4640 
    4641                 /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
    4642                 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
    4643         } //for (ig=0; ig<num_gauss; ig++
    4644 
    4645         /*Add Ke_gg to global matrix Kgg: */
    4646         MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
     4614                                        &Ke_g[0][0],0);
     4615
     4616                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke[i][j]+=Ke_g[i][j];
     4617        }
     4618
     4619        /*Add Ke to global matrix Kgg: */
     4620        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
    46474621               
    4648         xfree((void**)&first_gauss_area_coord);
    4649         xfree((void**)&second_gauss_area_coord);
    4650         xfree((void**)&third_gauss_area_coord);
    4651         xfree((void**)&gauss_weights);
     4622        /*Clean up*/
     4623        delete gauss;
    46524624        xfree((void**)&doflist);
    46534625}

  • issm/trunk/src/c/objects/Elements/TriaRef.cpp

    r4921 r5625  
    180180}
    181181/*}}}*/
    182 /*FUNCTION TriaRef::GetL {{{1*/
     182/*FUNCTION TriaRef::GetL(double* L, double* xyz_list, double* gauss,int numdof) {{{1*/
    183183void TriaRef::GetL(double* L, double* xyz_list, double* gauss,int numdof){
    184184        /*Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof.
    … …  
    221221}
    222222/*}}}*/
    223 /*FUNCTION TriaRef::GetJacobian {{{1*/
     223/*FUNCTION TriaRef::GetL(double* L, double* xyz_list,GaussTria* gauss,int numdof) {{{1*/
     224void TriaRef::GetL(double* L, double* xyz_list,GaussTria* gauss,int numdof){
     225        /*Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof.
     226         * For grid i, Li can be expressed in the actual coordinate system
     227         * by:
     228         *       numdof=1:
     229         *                 Li=h;
     230         *       numdof=2:
     231         *                 Li=[ h   0 ]
     232         *                    [ 0   h ]
     233         * where h is the interpolation function for grid i.
     234         *
     235         * We assume L has been allocated already, of size: numgrids (numdof=1), or numdofx(numdof*numgrids) (numdof=2)
     236         */
     237
     238        int i;
     239        const int NDOF2=2;
     240        const int numgrids=3;
     241        double l1l2l3[3];
     242
     243        /*Get l1l2l3 in actual coordinate system: */
     244        GetNodalFunctions(l1l2l3,gauss);
     245
     246        /*Build L: */
     247        if(numdof==1){
     248                for (i=0;i<numgrids;i++){
     249                        L[i]=l1l2l3[i];
     250                }
     251        }
     252        else{
     253                for (i=0;i<numgrids;i++){
     254                        *(L+numdof*numgrids*0+numdof*i)=l1l2l3[i]; //L[0][NDOF2*i]=dh1dh3[0][i];
     255                        *(L+numdof*numgrids*0+numdof*i+1)=0;
     256                        *(L+numdof*numgrids*1+numdof*i)=0;
     257                        *(L+numdof*numgrids*1+numdof*i+1)=l1l2l3[i];
     258                }
     259        }
     260}
     261/*}}}*/
     262/*FUNCTION TriaRef::GetJacobian(double* J, double* xyz_list,double* gauss) {{{1*/
    224263void TriaRef::GetJacobian(double* J, double* xyz_list,double* gauss){
    225264        /*The Jacobian is constant over the element, discard the gaussian points.
    … …  
    244283}
    245284/*}}}*/
    246 /*FUNCTION TriaRef::GetJacobianDeterminant2d {{{1*/
     285/*FUNCTION TriaRef::GetJacobian(double* J, double* xyz_list,GaussTria* gauss) {{{1*/
     286void TriaRef::GetJacobian(double* J, double* xyz_list,GaussTria* gauss){
     287        /*The Jacobian is constant over the element, discard the gaussian points.
     288         * J is assumed to have been allocated of size NDOF2xNDOF2.*/
     289
     290        const int NDOF2=2;
     291        const int numgrids=3;
     292        double x1,y1,x2,y2,x3,y3;
     293
     294        x1=*(xyz_list+numgrids*0+0);
     295        y1=*(xyz_list+numgrids*0+1);
     296        x2=*(xyz_list+numgrids*1+0);
     297        y2=*(xyz_list+numgrids*1+1);
     298        x3=*(xyz_list+numgrids*2+0);
     299        y3=*(xyz_list+numgrids*2+1);
     300
     301
     302        *(J+NDOF2*0+0)=0.5*(x2-x1);
     303        *(J+NDOF2*1+0)=SQRT3/6.0*(2*x3-x1-x2);
     304        *(J+NDOF2*0+1)=0.5*(y2-y1);
     305        *(J+NDOF2*1+1)=SQRT3/6.0*(2*y3-y1-y2);
     306}
     307/*}}}*/
     308/*FUNCTION TriaRef::GetJacobianDeterminant2d(double* Jdet, double* xyz_list,double* gauss) {{{1*/
    247309void TriaRef::GetJacobianDeterminant2d(double* Jdet, double* xyz_list,double* gauss){
    248310
     311        /*The Jacobian determinant is constant over the element, discard the gaussian points.
     312         * J is assumed to have been allocated of size NDOF2xNDOF2.*/
     313        double J[2][2];
     314
     315        /*Get Jacobian*/
     316        GetJacobian(&J[0][0],xyz_list,gauss);
     317
     318        /*Get Determinant*/
     319        Matrix2x2Determinant(Jdet,&J[0][0]);
     320        if(*Jdet<0) ISSMERROR("negative jacobian determinant!");
     321
     322}
     323/*}}}*/
     324/*FUNCTION TriaRef::GetJacobianDeterminant2d(double* Jdet, double* xyz_list,GaussTria* gauss) {{{1*/
     325void TriaRef::GetJacobianDeterminant2d(double* Jdet, double* xyz_list,GaussTria* gauss){
    249326        /*The Jacobian determinant is constant over the element, discard the gaussian points.
    250327         * J is assumed to have been allocated of size NDOF2xNDOF2.*/
    … …  
    296373}
    297374/*}}}*/
    298 /*FUNCTION TriaRef::GetNodalFunctions {{{1*/
     375/*FUNCTION TriaRef::GetNodalFunctions(double* l1l2l3, double* gauss) {{{1*/
    299376void TriaRef::GetNodalFunctions(double* l1l2l3, double* gauss){
    300377        /*This routine returns the values of the nodal functions  at the gaussian point.*/
    … …  
    303380        l1l2l3[1]=gauss[1];
    304381        l1l2l3[2]=gauss[2];
     382
     383}
     384/*}}}*/
     385/*FUNCTION TriaRef::GetNodalFunctions(double* l1l2l3,GaussTria* gauss) {{{1*/
     386void TriaRef::GetNodalFunctions(double* l1l2l3,GaussTria* gauss){
     387        /*This routine returns the values of the nodal functions  at the gaussian point.*/
     388
     389        l1l2l3[0]=gauss->coord1;
     390        l1l2l3[1]=gauss->coord2;
     391        l1l2l3[2]=gauss->coord3;
    305392
    306393}

  • issm/trunk/src/c/objects/Elements/TriaRef.h

    r4921 r5625  
    77#ifndef _TRIAREF_H_
    88#define _TRIAREF_H_
     9
     10class GaussTria;
    911
    1012class TriaRef{
    … …  
    2729                void GetBprimePrognostic(double* Bprime_prog, double* xyz_list, double* gauss);
    2830                void GetBPrognostic(double* B_prog, double* xyz_list, double* gauss);
    29                 void GetL(double* L, double* xyz_list, double* gauss,int numdof);
     31                void GetL(double* L, double* xyz_list,double* gauss,int numdof);
     32                void GetL(double* L, double* xyz_list,GaussTria* gauss,int numdof);
    3033                void GetJacobian(double* J, double* xyz_list,double* gauss);
     34                void GetJacobian(double* J, double* xyz_list,GaussTria* gauss);
    3135                void GetJacobianDeterminant2d(double* Jdet, double* xyz_list,double* gauss);
     36                void GetJacobianDeterminant2d(double* Jdet, double* xyz_list,GaussTria* gauss);
    3237                void GetJacobianDeterminant3d(double* Jdet, double* xyz_list,double* gauss);
    3338                void GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss);
    3439                void GetNodalFunctions(double* l1l2l3, double* gauss);
     40                void GetNodalFunctions(double* l1l2l3,GaussTria* gauss);
    3541                void GetNodalFunctionsDerivatives(double* l1l2l3,double* xyz_list, double* gauss);
    3642                void GetNodalFunctionsDerivativesReference(double* dl1dl3,double* gauss);

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

    r5578 r5625  
    2222/*Constraints: */
    2323#include "./Constraints/Spc.h"
     24
     25/*Gauss*/
     26#include "./Gauss/GaussTria.h"
    2427
    2528/*Loads: */
Note: See TracChangeset for help on using the changeset viewer.