Changeset 5831


Ignore:
Timestamp:
09/15/10 16:30:49 (15 years ago)
Author:
Mathieu Morlighem
Message:

Added merging of 2 matrices, to be completed (almost there)

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

Legend:

Unmodified
Added
Removed
  • issm/trunk/src/c/objects/Elements/Tria.cpp

    r5822 r5831  
    27552755void  Tria::CreateKMatrixDiagnosticMacAyeal(Mat Kgg,Mat Kff, Mat Kfs){
    27562756       
    2757         /*Intermediaries*/
    2758         double        *Ke_gg_viscous  = NULL;
    2759         double        *Ke_gg_friction = NULL;
    2760         ElementMatrix *Ke             = NULL;
    2761        
    2762         /*compute stiffness matrices on the g-set, at the element level: */
    2763         Ke_gg_viscous =  CreateKMatrixDiagnosticMacAyealViscous();
    2764         Ke_gg_friction = CreateKMatrixDiagnosticMacAyealFriction();
    2765 
    2766         /*Add Ke_gg values to Ke element stifness matrix: */
    2767         Ke=this->NewElementMatrix(MacAyealApproximationEnum);
    2768         if(Ke_gg_viscous) Ke->AddValues(Ke_gg_viscous);
    2769         if(Ke_gg_friction)Ke->AddValues(Ke_gg_friction);
     2757        /*compute all stiffness matrices for this element*/
     2758        ElementMatrix* Ke1=CreateKMatrixDiagnosticMacAyealViscous();
     2759        ElementMatrix* Ke2=CreateKMatrixDiagnosticMacAyealFriction();
     2760        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
    27702761
    27712762        /*Add Ke element stiffness matrix to global stiffness matrix: */
     
    27732764
    27742765        /*Free ressources:*/
     2766        delete Ke1;
     2767        delete Ke2;
    27752768        delete Ke;
    2776         xfree((void**)&Ke_gg_viscous);
    2777         xfree((void**)&Ke_gg_friction);
    27782769}
    27792770/*}}}*/
    27802771/*FUNCTION Tria::CreateKMatrixDiagnosticMacAyealViscous{{{1*/
    2781 double* Tria::CreateKMatrixDiagnosticMacAyealViscous(void){
    2782 
    2783         /* local declarations */
    2784         int             i,j;
    2785        
    2786         /*output: */
    2787         double*  Ke_gg=NULL;
    2788 
    2789         /* node data: */
    2790         const int    numdof=2*NUMVERTICES;
    2791         double       xyz_list[NUMVERTICES][3];
    2792 
    2793         /* gaussian points: */
    2794         int     ig;
    2795         GaussTria *gauss=NULL;
    2796 
    2797         /* material data: */
    2798         double viscosity; //viscosity
    2799         double newviscosity; //viscosity
    2800         double oldviscosity; //viscosity
    2801 
    2802         /* strain rate: */
    2803         double epsilon[3]; /* epsilon=[exx,eyy,exy];*/
    2804         double oldepsilon[3]; /* oldepsilon=[exx,eyy,exy];*/
    2805 
    2806         /* matrices: */
    2807         double B[3][numdof];
    2808         double Bprime[3][numdof];
    2809         double D[3][3]={0.0};  // material matrix, simple scalar matrix.
    2810         double D_scalar;
    2811 
    2812         /*parameters: */
    2813         double viscosity_overshoot;
    2814 
    2815         /* local element matrices: */
    2816         double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
    2817 
    2818         double Jdet;
    2819 
    2820         /*input parameters for structural analysis (diagnostic): */
    2821         double  thickness;
    2822        
    2823 
    2824         /*retrieve some parameters: */
    2825         this->parameters->FindParam(&viscosity_overshoot,ViscosityOvershootEnum);
    2826 
    2827         /*First, if we are on water, return empty matrix: */
    2828         if(IsOnWater()) return Ke_gg;
    2829 
    2830         /*allocate output: */
    2831         Ke_gg=(double*)xcalloc(numdof*numdof,sizeof(double));
    2832 
    2833         /* Get node coordinates and dof list: */
     2772ElementMatrix* Tria::CreateKMatrixDiagnosticMacAyealViscous(void){
     2773
     2774        /*Constants*/
     2775        const int  numdof=NDOF2*NUMVERTICES;
     2776
     2777        /*Intermediaries*/
     2778        int        i,j,ig;
     2779        double     xyz_list[NUMVERTICES][3];
     2780        GaussTria *gauss = NULL;
     2781        double     viscosity,newviscosity,oldviscosity;
     2782        double     viscosity_overshoot,thickness;
     2783        double     epsilon[3];    /* epsilon=[exx,eyy,exy];    */
     2784        double     oldepsilon[3]; /* oldepsilon=[exx,eyy,exy]; */
     2785        double     B[3][numdof];
     2786        double     Bprime[3][numdof];
     2787        double     D[3][3]   = {0.0};
     2788        double     D_scalar;
     2789        double     Ke_g[numdof][numdof];
     2790        double     Jdet;
     2791
     2792        /*Initialize Element matrix and return if necessary*/
     2793        if(IsOnWater()) return NULL;
     2794        ElementMatrix* Ke=this->NewElementMatrix(MacAyealApproximationEnum);
     2795
     2796        /*Retrieve all inputs and parameters*/
    28342797        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    2835 
    2836         /*Retrieve all inputs we will be needing: */
    28372798        Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input);
    28382799        Input* vx_input=inputs->GetInput(VxEnum);               ISSMASSERT(vx_input);
     
    28402801        Input* vxold_input=inputs->GetInput(VxOldEnum);         ISSMASSERT(vxold_input);
    28412802        Input* vyold_input=inputs->GetInput(VyOldEnum);         ISSMASSERT(vyold_input);
     2803        this->parameters->FindParam(&viscosity_overshoot,ViscosityOvershootEnum);
    28422804
    28432805        /* Start  looping on the number of gaussian points: */
     
    28472809                gauss->GaussPoint(ig);
    28482810
    2849                 /*Compute thickness at gaussian point: */
    28502811                thickness_input->GetParameterValue(&thickness, gauss);
    2851 
    2852                 /*Get strain rate from velocity: */
    28532812                this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    28542813                this->GetStrainRate2d(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
    2855 
    2856                 /*Get viscosity: */
    28572814                matice->GetViscosity2d(&viscosity, &epsilon[0]);
    28582815                matice->GetViscosity2d(&oldviscosity, &oldepsilon[0]);
    2859 
    2860                 /* Get Jacobian determinant: */
    28612816                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    28622817
    2863                 /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant
    2864                         onto this scalar matrix, so that we win some computational time: */
    28652818                newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
    28662819                D_scalar=2*newviscosity*thickness*gauss->weight*Jdet;
    2867 
    2868                 for (i=0;i<3;i++){
    2869                         D[i][i]=D_scalar;
    2870                 }
    2871 
    2872                 /*Get B and Bprime matrices: */
     2820                for (i=0;i<3;i++) D[i][i]=D_scalar;
     2821
    28732822                GetBMacAyeal(&B[0][0], &xyz_list[0][0], gauss);
    28742823                GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss);
    28752824
    2876                 /*  Do the triple product tB*D*Bprime: */
    2877                 TripleMultiply( &B[0][0],3,numdof,1,
     2825                TripleMultiply(&B[0][0],3,numdof,1,
    28782826                                        &D[0][0],3,3,0,
    28792827                                        &Bprime[0][0],3,numdof,0,
    2880                                         &Ke_gg_gaussian[0][0],0);
    2881 
    2882                 /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
    2883                 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) *(Ke_gg+numdof*j+i)+=Ke_gg_gaussian[i][j];
    2884 
     2828                                        &Ke_g[0][0],0);
     2829
     2830                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[numdof*j+i]+=Ke_g[i][j];
    28852831        }
    28862832
    28872833        /*Clean up and return*/
    28882834        delete gauss;
    2889         return Ke_gg;
     2835        return Ke;
    28902836}
    28912837/*}}}*/
    28922838/*FUNCTION Tria::CreateKMatrixDiagnosticMacAyealFriction {{{1*/
    28932839void  Tria::CreateKMatrixDiagnosticMacAyealFriction(Mat Kgg,Mat Kff, Mat Kfs){
    2894        
    2895         double        *Ke_gg_friction = NULL;
    2896         ElementMatrix *Ke             = NULL;
    2897        
    2898         /*compute stiffness matrices on the g-set, at the element level: */
    2899         Ke_gg_friction = CreateKMatrixDiagnosticMacAyealFriction();
    2900 
    2901         if(Ke_gg_friction){
    2902                 Ke=this->NewElementMatrix(MacAyealApproximationEnum);
    2903                 Ke->AddValues(Ke_gg_friction);
     2840
     2841        ElementMatrix* Ke=CreateKMatrixDiagnosticMacAyealFriction();
     2842
     2843        if(Ke){
    29042844                Ke->AddToGlobal(Kgg,Kff,Kfs);
    29052845                delete Ke;
    29062846        }
    2907 
    2908         /*Free ressources:*/
    2909         xfree((void**)&Ke_gg_friction);
    2910 }
    2911 
     2847}
    29122848/*}}}*/
    29132849/*FUNCTION Tria::CreateKMatrixDiagnosticMacAyealFriction {{{1*/
    2914 double*  Tria::CreateKMatrixDiagnosticMacAyealFriction(void){
    2915 
    2916         /* local declarations */
    2917         int       i,j;
    2918         int       analysis_type;
    2919 
    2920         /*output: */
    2921         double*   Ke_gg=NULL;
    2922 
    2923         /* node data: */
    2924         const int numdof   = 2 *NUMVERTICES;
    2925         double    xyz_list[NUMVERTICES][3];
    2926         int       numberofdofspernode=2;
    2927 
    2928         /* gaussian points: */
    2929         int     ig;
    2930         GaussTria *gauss=NULL;
    2931 
    2932         /* matrices: */
    2933         double L[2][numdof];
    2934         double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag
    2935         double DL_scalar;
    2936 
    2937         /* local element matrices: */
    2938         double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
    2939        
    2940         double Jdet;
    2941        
    2942         /*slope: */
    2943         double  slope[2]={0.0,0.0};
    2944         double  slope_magnitude;
    2945 
    2946         /*friction: */
    2947         Friction *friction = NULL;
    2948         double    alpha2;
    2949 
    2950         double MAXSLOPE=.06; // 6 %
    2951         double MOUNTAINKEXPONENT=10;
    2952 
    2953         /*inputs: */
    2954         int  drag_type;
    2955        
    2956 
    2957         /*retrive parameters: */
    2958         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    2959 
    2960         /*retrieve inputs :*/
    2961         inputs->GetParameterValue(&drag_type,DragTypeEnum);
     2850ElementMatrix* Tria::CreateKMatrixDiagnosticMacAyealFriction(void){
     2851
     2852        /*Constants*/
     2853        const int  numdof=NDOF2*NUMVERTICES;
     2854
     2855        /*Intermediaries*/
     2856        int        i,j,ig;
     2857        int        analysis_type,drag_type;
     2858        double     MAXSLOPE  = .06; // 6 %
     2859        double     MOUNTAINKEXPONENT = 10;
     2860        double     slope_magnitude,alpha2;
     2861        double     Jdet;
     2862        double     L[2][numdof];
     2863        double     DL[2][2]  = {{ 0,0 },{0,0}};
     2864        double     Ke_g[numdof][numdof];
     2865        double     DL_scalar;
     2866        double     slope[2]  = {0.0,0.0};
     2867        double     xyz_list[NUMVERTICES][3];
     2868        Friction  *friction = NULL;
     2869        GaussTria *gauss    = NULL;
     2870
     2871        /*Initialize Element matrix and return if necessary*/
     2872        if(IsOnWater() || IsOnShelf()) return NULL;
     2873        ElementMatrix* Ke=this->NewElementMatrix(MacAyealApproximationEnum);
     2874
     2875        /*Retrieve all inputs and parameters*/
     2876        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    29622877        Input* surface_input=inputs->GetInput(SurfaceEnum); ISSMASSERT(surface_input);
    29632878        Input* vx_input=inputs->GetInput(VxEnum);           ISSMASSERT(vx_input);
    29642879        Input* vy_input=inputs->GetInput(VyEnum);           ISSMASSERT(vy_input);
    29652880        Input* vz_input=inputs->GetInput(VzEnum);           ISSMASSERT(vz_input);
    2966        
    2967         /* Get node coordinates and dof list: */
    2968         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    2969 
    2970         if (IsOnShelf()){
    2971                 /*no friction, do nothing*/
    2972                 return Ke_gg;
    2973         }
    2974 
    2975         /*allocte Ke_gg matrix: */
    2976         Ke_gg=(double*)xcalloc(numdof*numdof,sizeof(double));
     2881        inputs->GetParameterValue(&drag_type,DragTypeEnum);
     2882        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    29772883
    29782884        /*build friction object, used later on: */
    2979         if (drag_type!=2)ISSMERROR(" non-viscous friction not supported yet!");
     2885        if (drag_type!=2) ISSMERROR(" non-viscous friction not supported yet!");
    29802886        friction=new Friction("2d",inputs,matpar,analysis_type);
    29812887
     
    29852891
    29862892                gauss->GaussPoint(ig);
    2987 
    2988                 /*Friction: */
    2989                 friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);
    29902893
    29912894                // If we have a slope > 6% for this element,  it means  we are on a mountain. In this particular case,
     
    29932896                surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
    29942897                slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));
    2995 
    2996                 if (slope_magnitude>MAXSLOPE){
    2997                         alpha2=pow((double)10,MOUNTAINKEXPONENT);
    2998                 }
    2999 
    3000                 /* Get Jacobian determinant: */
     2898                if(slope_magnitude>MAXSLOPE) alpha2=pow((double)10,MOUNTAINKEXPONENT);
     2899                else friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);
     2900
     2901                GetL(&L[0][0], &xyz_list[0][0], gauss,NDOF2);
    30012902                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    3002 
    3003                 /*Get L matrix: */
    3004                 GetL(&L[0][0], &xyz_list[0][0], gauss,numberofdofspernode);
    3005 
     2903                DL_scalar=alpha2*gauss->weight*Jdet;
     2904                for (i=0;i<2;i++) DL[i][i]=DL_scalar;
    30062905               
    3007                 DL_scalar=alpha2*gauss->weight*Jdet;
    3008                 for (i=0;i<2;i++){
    3009                         DL[i][i]=DL_scalar;
    3010                 }
    3011                
    3012                 /*  Do the triple producte tL*D*L: */
    30132906                TripleMultiply( &L[0][0],2,numdof,1,
    30142907                                        &DL[0][0],2,2,0,
    30152908                                        &L[0][0],2,numdof,0,
    3016                                         &Ke_gg_gaussian[0][0],0);
    3017 
    3018                 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) *(Ke_gg+numdof*i+j)+=Ke_gg_gaussian[i][j];
    3019 
     2909                                        &Ke_g[0][0],0);
     2910
     2911                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[numdof*j+i]+=Ke_g[i][j];
    30202912        }
    30212913
     
    30232915        delete gauss;
    30242916        delete friction;
    3025         return Ke_gg;
     2917        return Ke;
    30262918}
    30272919/*}}}*/
     
    59545846        int fsize;
    59555847        int ssize;
    5956         int* gexternaldoflist=NULL;
    5957         int* finternaldoflist=NULL;
    5958         int* fexternaldoflist=NULL;
    5959         int* sinternaldoflist=NULL;
    5960         int* sexternaldoflist=NULL;
     5848        int* gglobaldoflist=NULL;
     5849        int* flocaldoflist=NULL;
     5850        int* fglobaldoflist=NULL;
     5851        int* slocaldoflist=NULL;
     5852        int* sglobaldoflist=NULL;
    59615853        bool square=true;
    59625854
     
    59725864
    59735865        /*get dof lists for f and s set: */
    5974         if(!kff) gexternaldoflist=this->GetExternalDofList(approximation,GsetEnum);
    5975         else{
    5976                 finternaldoflist=this->GetInternalDofList(approximation,FsetEnum);
    5977                 fexternaldoflist=this->GetExternalDofList(approximation,FsetEnum);
    5978                 sinternaldoflist=this->GetInternalDofList(approximation,SsetEnum);
    5979                 sexternaldoflist=this->GetExternalDofList(approximation,SsetEnum);
     5866        gglobaldoflist=this->GetExternalDofList(approximation,GsetEnum);
     5867        if(kff){
     5868                flocaldoflist=this->GetInternalDofList(approximation,FsetEnum);
     5869                fglobaldoflist=this->GetExternalDofList(approximation,FsetEnum);
     5870                slocaldoflist=this->GetInternalDofList(approximation,SsetEnum);
     5871                sglobaldoflist=this->GetExternalDofList(approximation,SsetEnum);
    59805872        }
    59815873
    59825874        /*Use square constructor for ElementMatrix: */
    5983         if(!kff) Ke=new ElementMatrix(gsize,square,gexternaldoflist);
    5984         else     Ke=new ElementMatrix(gsize,square,finternaldoflist,fexternaldoflist,fsize,sinternaldoflist,sexternaldoflist,ssize);
     5875        if(!kff) Ke=new ElementMatrix(square,gglobaldoflist,gsize);
     5876        else     Ke=new ElementMatrix(square,gglobaldoflist,gsize,flocaldoflist,fglobaldoflist,fsize,slocaldoflist,sglobaldoflist,ssize);
    59855877
    59865878        /*Free ressources and return:*/
    5987         xfree((void**)&gexternaldoflist);
    5988         xfree((void**)&finternaldoflist);
    5989         xfree((void**)&fexternaldoflist);
    5990         xfree((void**)&sinternaldoflist);
    5991         xfree((void**)&sexternaldoflist);
     5879        xfree((void**)&gglobaldoflist);
     5880        xfree((void**)&flocaldoflist);
     5881        xfree((void**)&fglobaldoflist);
     5882        xfree((void**)&slocaldoflist);
     5883        xfree((void**)&sglobaldoflist);
    59925884
    59935885        return Ke;
     
    60045896        int gsize;
    60055897        int fsize;
    6006         int* gexternaldoflist=NULL;
    6007         int* finternaldoflist=NULL;
    6008         int* fexternaldoflist=NULL;
     5898        int* gglobaldoflist=NULL;
     5899        int* flocaldoflist=NULL;
     5900        int* fglobaldoflist=NULL;
    60095901
    60105902        /*retrieve some parameters: */
     
    60175909        /*get dof lists for f and s set: */
    60185910        if(!kff){
    6019                 gexternaldoflist=this->GetExternalDofList(approximation,GsetEnum);
     5911                gglobaldoflist=this->GetExternalDofList(approximation,GsetEnum);
    60205912        }
    60215913        else{
    6022                 finternaldoflist=this->GetInternalDofList(approximation,FsetEnum);
    6023                 fexternaldoflist=this->GetExternalDofList(approximation,FsetEnum);
     5914                flocaldoflist=this->GetInternalDofList(approximation,FsetEnum);
     5915                fglobaldoflist=this->GetExternalDofList(approximation,FsetEnum);
    60245916        }
    60255917
    60265918        /*constructor for ElementVector: */
    6027         if(!kff)pe=new ElementVector(gsize,gexternaldoflist);
    6028         else    pe=new ElementVector(gsize,finternaldoflist,fexternaldoflist,fsize);
     5919        if(!kff)pe=new ElementVector(gsize,gglobaldoflist);
     5920        else    pe=new ElementVector(gsize,flocaldoflist,fglobaldoflist,fsize);
    60295921
    60305922        /*Free ressources and return:*/
    6031         xfree((void**)&gexternaldoflist);
    6032         xfree((void**)&finternaldoflist);
    6033         xfree((void**)&fexternaldoflist);
     5923        xfree((void**)&gglobaldoflist);
     5924        xfree((void**)&flocaldoflist);
     5925        xfree((void**)&fglobaldoflist);
    60345926       
    60355927        return pe;
  • issm/trunk/src/c/objects/Elements/Tria.h

    r5780 r5831  
    125125                void      CreateKMatrixBalancedvelocities(Mat Kgg);
    126126                void      CreateKMatrixDiagnosticMacAyeal(Mat Kgg,Mat Kff, Mat Kfs);
    127                 double*  CreateKMatrixDiagnosticMacAyealViscous(void);
    128                 double*  CreateKMatrixDiagnosticMacAyealFriction(void);
    129                 void      CreateKMatrixDiagnosticMacAyealFriction(Mat Kgg,Mat Kff, Mat Kfs);
     127                ElementMatrix* CreateKMatrixDiagnosticMacAyealViscous(void);
     128                ElementMatrix* CreateKMatrixDiagnosticMacAyealFriction(void);
     129                void    CreateKMatrixDiagnosticMacAyealFriction(Mat Kgg,Mat Kff, Mat Kfs);
    130130                void      CreateKMatrixCouplingMacAyealPattynFriction(Mat Kgg);
    131131                void      CreateKMatrixDiagnosticPattynFriction(Mat Kgg);
  • issm/trunk/src/c/objects/Numerics/ElementMatrix.cpp

    r5827 r5831  
    4545}
    4646/*}}}*/
    47 /*FUNCTION ElementMatrix::ElementMatrix(int gsize,bool square,int* in_gglobaldoflist){{{1*/
    48 ElementMatrix::ElementMatrix(int gsize,bool square,int* in_gglobaldoflist){
     47/*FUNCTION ElementMatrix::ElementMatrix(ElementMatrix* Ke1, ElementMatrix* Ke2){{{1*/
     48ElementMatrix::ElementMatrix(ElementMatrix* Ke1, ElementMatrix* Ke2){
     49
     50        /*intermediaries*/
     51        int i,j,counter;
     52        int gsize,fsize,ssize;
     53        int* P=NULL;
     54        bool found;
     55
     56        /*If one of the two matrix is NULL, we copy the other one*/
     57        if(!Ke1 && !Ke2){
     58                ISSMERROR("Two input element matrices are NULL");
     59        }
     60        else if(!Ke1){
     61                this->Init(Ke2);
     62                return;
     63        }
     64        else if(!Ke2){
     65                this->Init(Ke1);
     66                return;
     67        }
     68
     69        /*General Case: Ke1 and Ke2 are not empty*/
     70        if(!Ke1->square || !Ke2->square) ISSMERROR("merging 2 non square matrices not implemented yet");
     71
     72        /*Initialize itransformation matrix Ke[P[i]] = Ke2[i]*/
     73        P=(int*)xmalloc(Ke2->nrows*sizeof(int));
     74
     75        /*1: Get the new numbering of Ke2 and get size of the new matrix*/
     76        gsize=Ke1->nrows;
     77        for(i=0;i<Ke2->nrows;i++){
     78                found=false;
     79                for(j=0;j<Ke1->nrows;j++){
     80                        if(Ke2->gglobaldoflist[i]==Ke1->gglobaldoflist[j]){
     81                                found=true; P[i]=j; break;
     82                        }
     83                }
     84                if(!found){
     85                        P[i]=gsize; gsize++;
     86                }
     87        }
     88
     89        /*2: Initialize static fields*/
     90        this->nrows=gsize;
     91        this->ncols=gsize;
     92        this->square=true;
     93        this->kff=Ke1->kff;
     94
     95        /*Gset and values*/
     96        this->gglobaldoflist=(int*)xmalloc(this->nrows*sizeof(int));
     97        this->values=(double*)xcalloc(this->nrows*this->ncols,sizeof(double));
     98        for(i=0;i<Ke1->nrows;i++){
     99                for(j=0;j<Ke1->ncols;j++){
     100                        this->values[i*this->ncols+j] += Ke1->values[i*Ke1->ncols+j];
     101                }
     102                this->gglobaldoflist[i]=Ke1->gglobaldoflist[i];
     103        }
     104        for(i=0;i<Ke2->nrows;i++){
     105                for(j=0;j<Ke2->ncols;j++){
     106                        this->values[P[i]*this->ncols+P[j]] += Ke2->values[i*Ke2->ncols+j];
     107                }
     108                this->gglobaldoflist[P[i]]=Ke2->gglobaldoflist[i];
     109        }
     110
     111        /*Fset*/
     112        fsize=Ke1->row_fsize;
     113        for(i=0;i<Ke2->row_fsize;i++){
     114                if(P[Ke2->row_flocaldoflist[i]] >= Ke1->nrows) fsize++;
     115        }
     116        printf("fsize = %i\n",fsize);
     117        ISSMERROR("STOP");
     118        if(fsize){
     119                this->row_flocaldoflist =(int*)xmalloc(fsize*sizeof(int));
     120                this->row_fglobaldoflist=(int*)xmalloc(fsize*sizeof(int));
     121                for(i=0;i<Ke1->row_fsize;i++){
     122                        this->row_flocaldoflist[i]=Ke1->row_flocaldoflist[i];
     123                }
     124                counter=Ke1->row_fsize;
     125                for(i=0;i<Ke2->row_fsize;i++){
     126                        if(P[Ke2->row_flocaldoflist[i]] >= Ke1->nrows){
     127                                this->row_flocaldoflist[counter]=P[Ke2->row_flocaldoflist[i]];
     128                        }
     129                }
     130        }
     131        else{
     132                this->row_flocaldoflist=NULL;
     133                this->row_fglobaldoflist=NULL;
     134        }
     135
     136        ssize=Ke1->row_ssize+Ke2->row_ssize;
     137        this->row_ssize=ssize;
     138        if(ssize){
     139                ISSMERROR("STOP");
     140                this->row_slocaldoflist=(int*)xmalloc(ssize*sizeof(int));
     141                this->row_sglobaldoflist=(int*)xmalloc(ssize*sizeof(int));
     142                //memcpy(this->row_slocaldoflist,slocaldoflist,ssize*sizeof(int));
     143                //memcpy(this->row_sglobaldoflist,sglobaldoflist,ssize*sizeof(int));
     144        }
     145        else{
     146                this->row_slocaldoflist=NULL;
     147                this->row_sglobaldoflist=NULL;
     148        }
     149
     150        /*don't do cols, we can pick them up from the rows: */
     151        this->col_fsize=0;
     152        this->col_flocaldoflist=NULL;
     153        this->col_fglobaldoflist=NULL;
     154        this->col_ssize=0;
     155        this->col_slocaldoflist=NULL;
     156        this->col_sglobaldoflist=NULL;
     157
     158        /*clean-up*/
     159        xfree((void**)&P);
     160}
     161/*}}}*/
     162/*FUNCTION ElementMatrix::ElementMatrix(bool square,int* in_gglobaldoflist,int gsize){{{1*/
     163ElementMatrix::ElementMatrix(bool square,int* in_gglobaldoflist,int gsize){
    49164
    50165        if(square=false)ISSMERROR(" calling square constructor with false square flag!");
     
    83198}
    84199/*}}}*/
    85 /*FUNCTION ElementMatrix::ElementMatrix(int gsize,bool square,int* flocaldoflist,int* fglobaldoflist,int fsize,int* slocaldoflist,int* sglobaldoflist,int ssize){{{1*/
    86 ElementMatrix::ElementMatrix(int gsize,bool square,int* flocaldoflist,int* fglobaldoflist,int fsize,int* slocaldoflist,int* sglobaldoflist,int ssize){
     200/*FUNCTION ElementMatrix::ElementMatrix(bool square,int* in_gglobaldoflist,int gsize,int* flocaldoflist,int* fglobaldoflist,int fsize,int* slocaldoflist,int* sglobaldoflist,int ssize){{{1*/
     201ElementMatrix::ElementMatrix(bool square,int* in_gglobaldoflist,int gsize,int* flocaldoflist,int* fglobaldoflist,int fsize,int* slocaldoflist,int* sglobaldoflist,int ssize){
    87202
    88203        if(square=false)ISSMERROR(" calling square constructor with false square flag!");
     204        ISSMASSERT(gsize);
    89205
    90206        this->nrows=gsize;
     
    95211        /*fill values with 0: */
    96212        this->values=(double*)xcalloc(this->nrows*this->ncols,sizeof(double));
    97        
     213
     214        /*g dofs*/
     215        this->gglobaldoflist=(int*)xmalloc(gsize*sizeof(int));
     216        memcpy(this->gglobaldoflist,in_gglobaldoflist,gsize*sizeof(int));
     217
    98218        /*row dofs: */
    99219        this->row_fsize=fsize;
     
    121241        }
    122242
    123         /*don't do cols, we can't pick them up from the rows: */
     243        /*don't do cols, we can pick them up from the rows: */
    124244        this->col_fsize=0;
    125245        this->col_flocaldoflist=NULL;
     
    128248        this->col_slocaldoflist=NULL;
    129249        this->col_sglobaldoflist=NULL;
    130 
    131         /*don't do g-set: */
    132         this->gglobaldoflist=NULL;
    133250
    134251}
     
    151268
    152269/*ElementMatrix specific routines: */
    153 /*FUNCTION ElementMatrix::AddValues(double* Ke_gg){{{1*/
     270/*FUNCTION ElementMatrix::AddValues{{{1*/
    154271void ElementMatrix::AddValues(double* Ke_gg){
    155272
     
    163280}
    164281/*}}}*/
    165 /*FUNCTION ElementMatrix::AddToGlobal(Mat Kgg, Mat Kff, Mat Kfs){{{1*/
     282/*FUNCTION ElementMatrix::AddToGlobal{{{1*/
    166283void ElementMatrix::AddToGlobal(Mat Kgg, Mat Kff, Mat Kfs){
    167284
     
    215332}
    216333/*}}}*/
    217 /*FUNCTION ElementMatrix::Echo(void){{{1*/
     334/*FUNCTION ElementMatrix::Echo{{{1*/
    218335void ElementMatrix::Echo(void){
    219336
     
    262379}
    263380/*}}}*/
     381/*FUNCTION ElementMatrix::Init{{{1*/
     382void ElementMatrix::Init(ElementMatrix* Ke){
     383
     384        ISSMASSERT(Ke);
     385
     386        this->nrows =Ke->nrows;
     387        this->ncols =Ke->ncols;
     388        this->square=Ke->square;
     389        this->kff   =Ke->kff;
     390        this->values=(double*)xmalloc(this->nrows*this->ncols*sizeof(double));
     391        memcpy(this->values,Ke->values,this->nrows*this->ncols*sizeof(double));
     392
     393        this->row_fsize=Ke->row_fsize;
     394        if(this->row_fsize){
     395                this->row_flocaldoflist=(int*)xmalloc(this->row_fsize*sizeof(int));
     396                memcpy(this->row_flocaldoflist,Ke->row_flocaldoflist,this->row_fsize*sizeof(int));
     397                this->row_fglobaldoflist=(int*)xmalloc(this->row_fsize*sizeof(int));
     398                memcpy(this->row_fglobaldoflist,Ke->row_fglobaldoflist,this->row_fsize*sizeof(int));
     399        }
     400        else{
     401                this->row_flocaldoflist=NULL;
     402                this->row_fglobaldoflist=NULL;
     403        }
     404
     405        this->row_ssize=Ke->row_ssize;
     406        if(this->row_ssize){
     407                this->row_slocaldoflist=(int*)xmalloc(this->row_ssize*sizeof(int));
     408                memcpy(this->row_slocaldoflist,Ke->row_slocaldoflist,this->row_ssize*sizeof(int));
     409                this->row_sglobaldoflist=(int*)xmalloc(this->row_ssize*sizeof(int));
     410                memcpy(this->row_sglobaldoflist,Ke->row_sglobaldoflist,this->row_ssize*sizeof(int));
     411        }
     412        else{
     413                this->row_slocaldoflist=NULL;
     414                this->row_sglobaldoflist=NULL;
     415        }
     416
     417        this->col_fsize=Ke->col_fsize;
     418        if(this->col_fsize){
     419                this->col_flocaldoflist=(int*)xmalloc(this->col_fsize*sizeof(int));
     420                memcpy(this->col_flocaldoflist,Ke->col_flocaldoflist,this->col_fsize*sizeof(int));
     421                this->col_fglobaldoflist=(int*)xmalloc(this->col_fsize*sizeof(int));
     422                memcpy(this->col_fglobaldoflist,Ke->col_fglobaldoflist,this->col_fsize*sizeof(int));
     423        }
     424        else{
     425                this->col_flocaldoflist=NULL;
     426                this->col_fglobaldoflist=NULL;
     427        }
     428
     429        this->col_ssize=Ke->col_ssize;
     430        if(this->col_ssize){
     431                this->col_slocaldoflist=(int*)xmalloc(this->col_ssize*sizeof(int));
     432                memcpy(this->col_slocaldoflist,Ke->col_slocaldoflist,this->col_ssize*sizeof(int));
     433                this->col_sglobaldoflist=(int*)xmalloc(this->col_ssize*sizeof(int));
     434                memcpy(this->col_sglobaldoflist,Ke->col_sglobaldoflist,this->col_ssize*sizeof(int));
     435        }
     436        else{
     437                this->col_slocaldoflist=NULL;
     438                this->col_sglobaldoflist=NULL;
     439        }
     440}
     441/*}}}*/
  • issm/trunk/src/c/objects/Numerics/ElementMatrix.h

    r5827 r5831  
    2121                int      nrows;
    2222                int      ncols;
    23                 double*  values;
    2423                bool     square;
    2524                bool     kff;
     25                double*  values;
    2626
    2727                //gset
     
    5050                /*ElementMatrix constructors, destructors {{{1*/
    5151                ElementMatrix();
    52                 ElementMatrix(int gsize,bool square,int* gglobaldoflist);
    53                 ElementMatrix(int gsize,bool square,int* flocaldoflist,int* fglobaldoflist,int fsize,int* slocaldoflist,int* sglobaldoflist,int ssize);
     52                ElementMatrix(ElementMatrix* Ke1,ElementMatrix* Ke2);
     53                ElementMatrix(bool square,int* gglobaldoflist,int gsize);
     54                ElementMatrix(bool square,int* gglobaldoflist,int gsize,int* flocaldoflist,int* fglobaldoflist,int fsize,int* slocaldoflist,int* sglobaldoflist,int ssize);
    5455                ~ElementMatrix();
    5556                /*}}}*/
     
    5859                void AddToGlobal(Mat Kgg, Mat Kff, Mat Kfs);
    5960                void Echo(void);
     61                void Init(ElementMatrix* Ke);
    6062                /*}}}*/
    6163};
Note: See TracChangeset for help on using the changeset viewer.