![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 5831
- Timestamp:
- 09/15/10 16:30:49 (15 years ago)
- Location:
- issm/trunk/src/c/objects
- Files:
-
- 4 edited
-
Elements/Tria.cpp (modified) (11 diffs)
-
Elements/Tria.h (modified) (1 diff)
-
Numerics/ElementMatrix.cpp (modified) (9 diffs)
-
Numerics/ElementMatrix.h (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Elements/Tria.cpp
r5822 r5831 2755 2755 void Tria::CreateKMatrixDiagnosticMacAyeal(Mat Kgg,Mat Kff, Mat Kfs){ 2756 2756 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); 2770 2761 2771 2762 /*Add Ke element stiffness matrix to global stiffness matrix: */ … … 2773 2764 2774 2765 /*Free ressources:*/ 2766 delete Ke1; 2767 delete Ke2; 2775 2768 delete Ke; 2776 xfree((void**)&Ke_gg_viscous);2777 xfree((void**)&Ke_gg_friction);2778 2769 } 2779 2770 /*}}}*/ 2780 2771 /*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: */ 2772 ElementMatrix* 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*/ 2834 2797 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 2835 2836 /*Retrieve all inputs we will be needing: */2837 2798 Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); 2838 2799 Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); … … 2840 2801 Input* vxold_input=inputs->GetInput(VxOldEnum); ISSMASSERT(vxold_input); 2841 2802 Input* vyold_input=inputs->GetInput(VyOldEnum); ISSMASSERT(vyold_input); 2803 this->parameters->FindParam(&viscosity_overshoot,ViscosityOvershootEnum); 2842 2804 2843 2805 /* Start looping on the number of gaussian points: */ … … 2847 2809 gauss->GaussPoint(ig); 2848 2810 2849 /*Compute thickness at gaussian point: */2850 2811 thickness_input->GetParameterValue(&thickness, gauss); 2851 2852 /*Get strain rate from velocity: */2853 2812 this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); 2854 2813 this->GetStrainRate2d(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input); 2855 2856 /*Get viscosity: */2857 2814 matice->GetViscosity2d(&viscosity, &epsilon[0]); 2858 2815 matice->GetViscosity2d(&oldviscosity, &oldepsilon[0]); 2859 2860 /* Get Jacobian determinant: */2861 2816 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); 2862 2817 2863 /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant2864 onto this scalar matrix, so that we win some computational time: */2865 2818 newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); 2866 2819 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 2873 2822 GetBMacAyeal(&B[0][0], &xyz_list[0][0], gauss); 2874 2823 GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss); 2875 2824 2876 /* Do the triple product tB*D*Bprime: */ 2877 TripleMultiply( &B[0][0],3,numdof,1, 2825 TripleMultiply(&B[0][0],3,numdof,1, 2878 2826 &D[0][0],3,3,0, 2879 2827 &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]; 2885 2831 } 2886 2832 2887 2833 /*Clean up and return*/ 2888 2834 delete gauss; 2889 return Ke _gg;2835 return Ke; 2890 2836 } 2891 2837 /*}}}*/ 2892 2838 /*FUNCTION Tria::CreateKMatrixDiagnosticMacAyealFriction {{{1*/ 2893 2839 void 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){ 2904 2844 Ke->AddToGlobal(Kgg,Kff,Kfs); 2905 2845 delete Ke; 2906 2846 } 2907 2908 /*Free ressources:*/ 2909 xfree((void**)&Ke_gg_friction); 2910 } 2911 2847 } 2912 2848 /*}}}*/ 2913 2849 /*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); 2850 ElementMatrix* 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); 2962 2877 Input* surface_input=inputs->GetInput(SurfaceEnum); ISSMASSERT(surface_input); 2963 2878 Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); 2964 2879 Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); 2965 2880 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); 2977 2883 2978 2884 /*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!"); 2980 2886 friction=new Friction("2d",inputs,matpar,analysis_type); 2981 2887 … … 2985 2891 2986 2892 gauss->GaussPoint(ig); 2987 2988 /*Friction: */2989 friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);2990 2893 2991 2894 // If we have a slope > 6% for this element, it means we are on a mountain. In this particular case, … … 2993 2896 surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss); 2994 2897 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); 3001 2902 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; 3006 2905 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: */3013 2906 TripleMultiply( &L[0][0],2,numdof,1, 3014 2907 &DL[0][0],2,2,0, 3015 2908 &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]; 3020 2912 } 3021 2913 … … 3023 2915 delete gauss; 3024 2916 delete friction; 3025 return Ke _gg;2917 return Ke; 3026 2918 } 3027 2919 /*}}}*/ … … 5954 5846 int fsize; 5955 5847 int ssize; 5956 int* g externaldoflist=NULL;5957 int* f internaldoflist=NULL;5958 int* f externaldoflist=NULL;5959 int* s internaldoflist=NULL;5960 int* s externaldoflist=NULL;5848 int* gglobaldoflist=NULL; 5849 int* flocaldoflist=NULL; 5850 int* fglobaldoflist=NULL; 5851 int* slocaldoflist=NULL; 5852 int* sglobaldoflist=NULL; 5961 5853 bool square=true; 5962 5854 … … 5972 5864 5973 5865 /*get dof lists for f and s set: */ 5974 if(!kff) gexternaldoflist=this->GetExternalDofList(approximation,GsetEnum);5975 else{5976 f internaldoflist=this->GetInternalDofList(approximation,FsetEnum);5977 f externaldoflist=this->GetExternalDofList(approximation,FsetEnum);5978 s internaldoflist=this->GetInternalDofList(approximation,SsetEnum);5979 s externaldoflist=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); 5980 5872 } 5981 5873 5982 5874 /*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); 5985 5877 5986 5878 /*Free ressources and return:*/ 5987 xfree((void**)&g externaldoflist);5988 xfree((void**)&f internaldoflist);5989 xfree((void**)&f externaldoflist);5990 xfree((void**)&s internaldoflist);5991 xfree((void**)&s externaldoflist);5879 xfree((void**)&gglobaldoflist); 5880 xfree((void**)&flocaldoflist); 5881 xfree((void**)&fglobaldoflist); 5882 xfree((void**)&slocaldoflist); 5883 xfree((void**)&sglobaldoflist); 5992 5884 5993 5885 return Ke; … … 6004 5896 int gsize; 6005 5897 int fsize; 6006 int* g externaldoflist=NULL;6007 int* f internaldoflist=NULL;6008 int* f externaldoflist=NULL;5898 int* gglobaldoflist=NULL; 5899 int* flocaldoflist=NULL; 5900 int* fglobaldoflist=NULL; 6009 5901 6010 5902 /*retrieve some parameters: */ … … 6017 5909 /*get dof lists for f and s set: */ 6018 5910 if(!kff){ 6019 g externaldoflist=this->GetExternalDofList(approximation,GsetEnum);5911 gglobaldoflist=this->GetExternalDofList(approximation,GsetEnum); 6020 5912 } 6021 5913 else{ 6022 f internaldoflist=this->GetInternalDofList(approximation,FsetEnum);6023 f externaldoflist=this->GetExternalDofList(approximation,FsetEnum);5914 flocaldoflist=this->GetInternalDofList(approximation,FsetEnum); 5915 fglobaldoflist=this->GetExternalDofList(approximation,FsetEnum); 6024 5916 } 6025 5917 6026 5918 /*constructor for ElementVector: */ 6027 if(!kff)pe=new ElementVector(gsize,g externaldoflist);6028 else pe=new ElementVector(gsize,f internaldoflist,fexternaldoflist,fsize);5919 if(!kff)pe=new ElementVector(gsize,gglobaldoflist); 5920 else pe=new ElementVector(gsize,flocaldoflist,fglobaldoflist,fsize); 6029 5921 6030 5922 /*Free ressources and return:*/ 6031 xfree((void**)&g externaldoflist);6032 xfree((void**)&f internaldoflist);6033 xfree((void**)&f externaldoflist);5923 xfree((void**)&gglobaldoflist); 5924 xfree((void**)&flocaldoflist); 5925 xfree((void**)&fglobaldoflist); 6034 5926 6035 5927 return pe; -
issm/trunk/src/c/objects/Elements/Tria.h
r5780 r5831 125 125 void CreateKMatrixBalancedvelocities(Mat Kgg); 126 126 void CreateKMatrixDiagnosticMacAyeal(Mat Kgg,Mat Kff, Mat Kfs); 127 double*CreateKMatrixDiagnosticMacAyealViscous(void);128 double*CreateKMatrixDiagnosticMacAyealFriction(void);129 void 127 ElementMatrix* CreateKMatrixDiagnosticMacAyealViscous(void); 128 ElementMatrix* CreateKMatrixDiagnosticMacAyealFriction(void); 129 void CreateKMatrixDiagnosticMacAyealFriction(Mat Kgg,Mat Kff, Mat Kfs); 130 130 void CreateKMatrixCouplingMacAyealPattynFriction(Mat Kgg); 131 131 void CreateKMatrixDiagnosticPattynFriction(Mat Kgg); -
issm/trunk/src/c/objects/Numerics/ElementMatrix.cpp
r5827 r5831 45 45 } 46 46 /*}}}*/ 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*/ 48 ElementMatrix::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*/ 163 ElementMatrix::ElementMatrix(bool square,int* in_gglobaldoflist,int gsize){ 49 164 50 165 if(square=false)ISSMERROR(" calling square constructor with false square flag!"); … … 83 198 } 84 199 /*}}}*/ 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*/ 201 ElementMatrix::ElementMatrix(bool square,int* in_gglobaldoflist,int gsize,int* flocaldoflist,int* fglobaldoflist,int fsize,int* slocaldoflist,int* sglobaldoflist,int ssize){ 87 202 88 203 if(square=false)ISSMERROR(" calling square constructor with false square flag!"); 204 ISSMASSERT(gsize); 89 205 90 206 this->nrows=gsize; … … 95 211 /*fill values with 0: */ 96 212 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 98 218 /*row dofs: */ 99 219 this->row_fsize=fsize; … … 121 241 } 122 242 123 /*don't do cols, we can 'tpick them up from the rows: */243 /*don't do cols, we can pick them up from the rows: */ 124 244 this->col_fsize=0; 125 245 this->col_flocaldoflist=NULL; … … 128 248 this->col_slocaldoflist=NULL; 129 249 this->col_sglobaldoflist=NULL; 130 131 /*don't do g-set: */132 this->gglobaldoflist=NULL;133 250 134 251 } … … 151 268 152 269 /*ElementMatrix specific routines: */ 153 /*FUNCTION ElementMatrix::AddValues (double* Ke_gg){{{1*/270 /*FUNCTION ElementMatrix::AddValues{{{1*/ 154 271 void ElementMatrix::AddValues(double* Ke_gg){ 155 272 … … 163 280 } 164 281 /*}}}*/ 165 /*FUNCTION ElementMatrix::AddToGlobal (Mat Kgg, Mat Kff, Mat Kfs){{{1*/282 /*FUNCTION ElementMatrix::AddToGlobal{{{1*/ 166 283 void ElementMatrix::AddToGlobal(Mat Kgg, Mat Kff, Mat Kfs){ 167 284 … … 215 332 } 216 333 /*}}}*/ 217 /*FUNCTION ElementMatrix::Echo (void){{{1*/334 /*FUNCTION ElementMatrix::Echo{{{1*/ 218 335 void ElementMatrix::Echo(void){ 219 336 … … 262 379 } 263 380 /*}}}*/ 381 /*FUNCTION ElementMatrix::Init{{{1*/ 382 void 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 21 21 int nrows; 22 22 int ncols; 23 double* values;24 23 bool square; 25 24 bool kff; 25 double* values; 26 26 27 27 //gset … … 50 50 /*ElementMatrix constructors, destructors {{{1*/ 51 51 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); 54 55 ~ElementMatrix(); 55 56 /*}}}*/ … … 58 59 void AddToGlobal(Mat Kgg, Mat Kff, Mat Kfs); 59 60 void Echo(void); 61 void Init(ElementMatrix* Ke); 60 62 /*}}}*/ 61 63 };
Note:
See TracChangeset
for help on using the changeset viewer.
