Changeset 16371

Timestamp:

10/10/13 15:58:56 (11 years ago)

Author:

Mathieu Morlighem

Message:

NEW: preparing mesh deformation module

Location:

issm/trunk-jpl

Files:

• m4/issm_options.m4 (modified) (1 diff)
• src/c/Makefile.am (modified) (2 diffs)
• src/c/analyses/analyses.h (modified) (1 diff)
• src/c/analyses/meshdeformation_core.cpp (added)
• src/c/classes/Elements/Tria.cpp (modified) (2 diffs)
• src/c/classes/Elements/TriaRef.cpp (modified) (17 diffs)
• src/c/classes/Elements/TriaRef.h (modified) (2 diffs)
• src/c/modules/ModelProcessorx/MeshDeformation (added)
• src/c/modules/ModelProcessorx/MeshDeformation/CreateConstraintsMeshDeformation.cpp (added)
• src/c/modules/ModelProcessorx/MeshDeformation/CreateLoadsMeshDeformation.cpp (added)
• src/c/modules/ModelProcessorx/MeshDeformation/CreateNodesMeshDeformation.cpp (added)
• src/c/modules/ModelProcessorx/MeshDeformation/UpdateElementsMeshDeformation.cpp (added)
• src/c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
• src/c/shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• src/c/shared/Enum/StringToEnumx.cpp (modified) (3 diffs)
• src/m/enum/EnumDefinitions.py (modified) (1 diff)
• src/m/enum/MeshdeformationAnalysisEnum.m (added)
• src/m/enum/MeshdeformationSolutionEnum.m (added)

issm/trunk-jpl/m4/issm_options.m4

@@ -1766 +1766 @@
         AC_MSG_RESULT($HAVE_SLOPE)
         dnl }}}
+        dnl with-meshdeformation{{{
+                AC_ARG_WITH([meshdeformation],
+                                        AS_HELP_STRING([--with-meshdeformation = YES], [compile with meshdeformation capabilities (default is yes)]),
+                                        [MESHDEFORMATION=$withval],[MESHDEFORMATION=yes]) 
+                  AC_MSG_CHECKING(for meshdeformation capability compilation)
+
+                  if test "x$MESHDEFORMATION" = "xyes"; then
+
+                        dnl defaults
+                          HAVE_MESHDEFORMATION=yes
+
+                          AC_DEFINE([_HAVE_MESHDEFORMATION_],[1],[with meshdeformation capability])
+                  else
+                                HAVE_MESHDEFORMATION=no
+                  fi
+          AM_CONDITIONAL([MESHDEFORMATION], [test x$HAVE_MESHDEFORMATION = xyes])
+          AC_MSG_RESULT($HAVE_MESHDEFORMATION)
+          dnl }}}
         dnl with-groundingline{{{
         AC_ARG_WITH([groundingline],

issm/trunk-jpl/src/c/Makefile.am

@@ -528 +528 @@
                                           ./analyses/bedslope_core.cpp
 #}}}
+#MeshDeformation sources  {{{
+meshdeformation_sources =  ./modules/ModelProcessorx/MeshDeformation/UpdateElementsMeshDeformation.cpp\
+                                          ./modules/ModelProcessorx/MeshDeformation/CreateNodesMeshDeformation.cpp \
+                                          ./modules/ModelProcessorx/MeshDeformation/CreateConstraintsMeshDeformation.cpp\
+                                          ./modules/ModelProcessorx/MeshDeformation/CreateLoadsMeshDeformation.cpp\
+                                          ./analyses/meshdeformation_core.cpp
+#}}}
 #Gia sources  {{{
 gia_sources =  ./analyses/gia_core.cpp\
@@ -934 +941 @@
 endif
 
+if MESHDEFORMATION
+issm_sources  +=  $(meshdeformation_sources)
+endif
+
 if GROUNDINGLINE
 issm_sources +=  $(groundingline_sources)

issm/trunk-jpl/src/c/analyses/analyses.h

@@ -26 +26 @@
 void surfaceslope_core(FemModel* femmodel);
 void bedslope_core(FemModel* femmodel);
+void meshdeformation_core(FemModel* femmodel);
 void control_core(FemModel* femmodel);
 void controltao_core(FemModel* femmodel);

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -3209 +3209 @@
         int         i,j;
         int         analysis_type,approximation;
-        IssmDouble  alpha2,Jdet2d;
+        int         indices[2];
+        IssmDouble  alpha2,Jdet;
         IssmDouble  FSreconditioning,viscosity;
         IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
@@ -3215 +3216 @@
         IssmDouble  xyz_list_seg[NUMVERTICES1D][3];
         Friction   *friction = NULL;
-        GaussTria  *gauss    = NULL;
 
         /*Initialize Element matrix and return if necessary*/
         if(IsFloating() || !HasEdgeOnBed()) return NULL;
 
-        _error_("STOP");
+        /*Fetch number of nodes and dof for this finite element*/
+        int vnumnodes = this->NumberofNodesVelocity();
+        int pnumnodes = this->NumberofNodesPressure();
+        int  numdof   = vnumnodes*NDOF2 + pnumnodes*NDOF1;
+
+        /*Prepare coordinate system list*/
+        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
+        for(i=0;i<vnumnodes;i++) cs_list[i]           = XYEnum;
+        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
+
+        /*Initialize Element matrix and vectors*/
+        ElementMatrix* Ke        = new ElementMatrix(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
+        IssmDouble*    BFriction = xNew<IssmDouble>(2*numdof);
+        IssmDouble*    D         = xNewZeroInit<IssmDouble>(2*2);
+
+        /*Retrieve all inputs and parameters*/
+        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+        parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
+        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
+        Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
+
+        /*Get vertex indices that lie on bed*/
+        this->EdgeOnBedIndices(&indices[0],&indices[1]);
+        for(i=0;i<NUMVERTICES1D;i++) for(j=0;j<2;j++) xyz_list_seg[i][j]=xyz_list[indices[i]][j];
+
+        /*build friction object, used later on: */
+        friction=new Friction("1d",inputs,matpar,analysis_type);
+
+        /* Start looping on the number of gauss 1d (nodes on the bedrock) */
+        GaussTria* gauss=new GaussTria(indices[0],indices[1],2);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+
+                GetSegmentJacobianDeterminant(&Jdet,&xyz_list_seg[0][0],gauss);
+
+                _error_("STOP");
+        }
+
+        /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/
+        //TransformStiffnessMatrixCoord(Ke,nodes,(vnumnodes+pnumnodes),cs_list);
+
+        /*Clean up and return*/
+        delete gauss;
+        delete friction;
+        xDelete<IssmDouble>(BFriction);
+        xDelete<IssmDouble>(D);
+        xDelete<int>(cs_list);
+        return Ke;
 }
 /*}}}*/

issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp

@@ -57 +57 @@
          *       Bi=[ dN/dx ]
          *          [ dN/dy ]
-         * where N is the interpolation function for node i.
+         * where N is the finiteelement function for node i.
          *
          * We assume B has been allocated already, of size: 3x(NDOF2*numnodes)
@@ -87 +87 @@
          *          [   0           dN/dy  ]
          *          [ 1/2*dN/dy  1/2*dN/dx ]
-         * where N is the interpolation function for node i.
+         * where N is the finiteelement function for node i.
          *
          * We assume B has been allocated already, of size: 3x(NDOF2*numnodes)
@@ -122 +122 @@
          *          [       0       dN/dy   ]
          *          [  1/2*dN/dy  1/2*dN/dx ]
-         * where N is the interpolation function for node i.
+         * where N is the finiteelement function for node i.
          *
          * We assume B has been allocated already, of size: 3x(NDOF2*numnodes)
@@ -215 +215 @@
          *                                      [ phi_p ]
          *                                      [  0    ]
-         *      where phi is the interpolation function for node i.
+         *      where phi is the finiteelement function for node i.
          *      Same thing for Bb except the last column that does not exist.
          */
@@ -271 +271 @@
          *                                      [  0  ]
          *                                      [ phi ]
-         *      where phi is the interpolation function for node i.
+         *      where phi is the finiteelement function for node i.
          */
 
@@ -317 +317 @@
          *       Bi=[ N ]
          *          [ N ]
-         * where N is the interpolation function for node i.
+         * where N is the finiteelement function for node i.
          *
          * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes)
@@ -348 +348 @@
          *                [   dN/dx  2*dN/dy ]
          *                [   dN/dy    dN/dx ]
-         * where hNis the interpolation function for node i.
+         * where hNis the finiteelement function for node i.
          *
          * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
@@ -383 +383 @@
          *               [  dN/dy  dN/dx ]
          N               [  dN/dx  dN/dy ]
-         * where N is the interpolation function for node i.
+         * where N is the finiteelement function for node i.
          *
          * We assume Bprime has been allocated already, of size: 3x(NDOF2*numnodes)
@@ -418 +418 @@
          *       Bi_prime=[ dN/dx ]
          *                [ dN/dy ]
-         * where N is the interpolation function for node i.
+         * where N is the finiteelement function for node i.
          *
          * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
@@ -447 +447 @@
          *                 Bi=[ N   0 ]
          *                    [ 0   N ]
-         * where N is the interpolation function for node i.
+         * where N is the finiteelement function for node i.
          *
          * We assume B has been allocated already, of size: 2 x (numdof*numnodes)
@@ -542 +542 @@
 }
 /*}}}*/
-/*FUNCTION TriaRef::GetNodalFunctions(IssmDouble* basis,GaussTria* gauss,int interpolation){{{*/
-void TriaRef::GetNodalFunctions(IssmDouble* basis,GaussTria* gauss,int interpolation){
+/*FUNCTION TriaRef::GetNodalFunctions(IssmDouble* basis,GaussTria* gauss,int finiteelement){{{*/
+void TriaRef::GetNodalFunctions(IssmDouble* basis,GaussTria* gauss,int finiteelement){
         /*This routine returns the values of the nodal functions  at the gaussian point.*/
 
         _assert_(basis);
 
-        switch(interpolation){
+        switch(finiteelement){
                 case P1Enum: case P1DGEnum:
                         basis[0]=gauss->coord1;
@@ -687 +687 @@
 }
 /*}}}*/
-/*FUNCTION TriaRef::GetNodalFunctionsDerivatives{{{*/
+/*FUNCTION TriaRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTria* gauss){{{*/
 void TriaRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTria* gauss){
+
+        GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss,this->element_type);
+
+}
+/*}}}*/
+/*FUNCTION TriaRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTria* gauss,int finiteelement){{{*/
+void TriaRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTria* gauss,int finiteelement){
 
         /*This routine returns the values of the nodal functions derivatives  (with respect to the 
@@ -695 +702 @@
 
         /*Fetch number of nodes for this finite element*/
-        int numnodes = this->NumberofNodes();
+        int numnodes = this->NumberofNodes(finiteelement);
 
         /*Get nodal functions derivatives in reference triangle*/
         IssmDouble* dbasis_ref=xNew<IssmDouble>(2*numnodes);
-        GetNodalFunctionsDerivativesReference(dbasis_ref,gauss); 
+        GetNodalFunctionsDerivativesReference(dbasis_ref,gauss,finiteelement); 
 
         /*Get Jacobian invert: */
@@ -784 +791 @@
 }
 /*}}}*/
-/*FUNCTION TriaRef::GetNodalFunctionsDerivativesReference{{{*/
+/*FUNCTION TriaRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTria* gauss){{{*/
 void TriaRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTria* gauss){
         /*This routine returns the values of the nodal functions derivatives  (with respect to the 
          * natural coordinate system) at the gaussian point. */
 
+        GetNodalFunctionsDerivativesReference(dbasis,gauss,this->element_type);
+
+}
+/*}}}*/
+/*FUNCTION TriaRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTria* gauss,int finiteelement){{{*/
+void TriaRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTria* gauss,int finiteelement){
+        /*This routine returns the values of the nodal functions derivatives  (with respect to the 
+         * natural coordinate system) at the gaussian point. */
+
         _assert_(dbasis && gauss);
 
-        switch(this->element_type){
+        switch(finiteelement){
                 case P1Enum: case P1DGEnum:
                         /*Nodal function 1*/
@@ -882 +898 @@
 }
 /*}}}*/
-/*FUNCTION TriaRef::GetInputValue(IssmDouble* p, IssmDouble* plist, GaussTria* gauss,int interpolation){{{*/
-void TriaRef::GetInputValue(IssmDouble* p, IssmDouble* plist, GaussTria* gauss,int interpolation){
+/*FUNCTION TriaRef::GetInputValue(IssmDouble* p, IssmDouble* plist, GaussTria* gauss,int finiteelement){{{*/
+void TriaRef::GetInputValue(IssmDouble* p, IssmDouble* plist, GaussTria* gauss,int finiteelement){
 
         /*Output*/
@@ -889 +905 @@
 
         /*Fetch number of nodes for this finite element*/
-        int numnodes = this->NumberofNodes(interpolation);
+        int numnodes = this->NumberofNodes(finiteelement);
 
         /*Get nodal functions*/
         IssmDouble* basis=xNew<IssmDouble>(numnodes);
-        GetNodalFunctions(basis, gauss,interpolation);
+        GetNodalFunctions(basis, gauss,finiteelement);
 
         /*Calculate parameter for this Gauss point*/
@@ -909 +925 @@
 }
 /*}}}*/
-/*FUNCTION TriaRef::NumberofNodes(int interpolation){{{*/
-int TriaRef::NumberofNodes(int interpolation){
-
-        switch(interpolation){
+/*FUNCTION TriaRef::NumberofNodes(int finiteelement){{{*/
+int TriaRef::NumberofNodes(int finiteelement){
+
+        switch(finiteelement){
                 case P1Enum:                return NUMNODESP1;
                 case P1DGEnum:              return NUMNODESP1;

issm/trunk-jpl/src/c/classes/Elements/TriaRef.h

@@ -38 +38 @@
                 void GetJacobianInvert(IssmDouble*  Jinv, IssmDouble* xyz_list,GaussTria* gauss);
                 void GetNodalFunctions(IssmDouble* basis,GaussTria* gauss);
-                void GetNodalFunctions(IssmDouble* basis,GaussTria* gauss,int interpolation);
+                void GetNodalFunctions(IssmDouble* basis,GaussTria* gauss,int finiteelement);
                 void GetNodalFunctionsVelocity(IssmDouble* basis, GaussTria* gauss);
                 void GetNodalFunctionsPressure(IssmDouble* basis, GaussTria* gauss);
@@ -45 +45 @@
                 void GetSegmentBprimeFlux(IssmDouble* Bprime,GaussTria* gauss, int index1,int index2);
                 void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTria* gauss);
+                void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTria* gauss,int finiteelement);
                 void GetNodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,GaussTria* gauss);
                 void GetNodalFunctionsDerivativesPressure(IssmDouble* dbasis,IssmDouble* xyz_list,GaussTria* gauss);
                 void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTria* gauss);
+                void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTria* gauss,int finiteelement);
                 void GetInputValue(IssmDouble* pp, IssmDouble* plist, GaussTria* gauss);
-                void GetInputValue(IssmDouble* pp, IssmDouble* plist, GaussTria* gauss,int interpolation);
+                void GetInputValue(IssmDouble* pp, IssmDouble* plist, GaussTria* gauss,int finiteelement);
                 void GetInputDerivativeValue(IssmDouble* pp, IssmDouble* plist,IssmDouble* xyz_list, GaussTria* gauss);
 
                 int  NumberofNodes(void);
-                int  NumberofNodes(int interpolation);
+                int  NumberofNodes(int finiteelement);
                 int  NumberofNodesVelocity(void);
                 int  NumberofNodesPressure(void);

issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

@@ -332 +332 @@
         GiaSolutionEnum,
         GiaAnalysisEnum,
+        MeshdeformationSolutionEnum,
+        MeshdeformationAnalysisEnum,
         /*}}}*/
         /*Approximations {{{*/

issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

@@ -338 +338 @@
                 case GiaSolutionEnum : return "GiaSolution";
                 case GiaAnalysisEnum : return "GiaAnalysis";
+                case MeshdeformationSolutionEnum : return "MeshdeformationSolution";
+                case MeshdeformationAnalysisEnum : return "MeshdeformationAnalysis";
                 case ApproximationEnum : return "Approximation";
                 case NoneApproximationEnum : return "NoneApproximation";

issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

@@ -344 +344 @@
               else if (strcmp(name,"GiaSolution")==0) return GiaSolutionEnum;
               else if (strcmp(name,"GiaAnalysis")==0) return GiaAnalysisEnum;
+              else if (strcmp(name,"MeshdeformationSolution")==0) return MeshdeformationSolutionEnum;
+              else if (strcmp(name,"MeshdeformationAnalysis")==0) return MeshdeformationAnalysisEnum;
               else if (strcmp(name,"Approximation")==0) return ApproximationEnum;
               else if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum;
@@ -381 +383 @@
               else if (strcmp(name,"ExternalResult")==0) return ExternalResultEnum;
               else if (strcmp(name,"FileParam")==0) return FileParamEnum;
-              else if (strcmp(name,"Input")==0) return InputEnum;
-              else if (strcmp(name,"IntInput")==0) return IntInputEnum;
          else stage=4;
    }
    if(stage==4){
-              if (strcmp(name,"IntParam")==0) return IntParamEnum;
+              if (strcmp(name,"Input")==0) return InputEnum;
+              else if (strcmp(name,"IntInput")==0) return IntInputEnum;
+              else if (strcmp(name,"IntParam")==0) return IntParamEnum;
               else if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum;
               else if (strcmp(name,"TransientParam")==0) return TransientParamEnum;
@@ -504 +506 @@
               else if (strcmp(name,"StressTensor")==0) return StressTensorEnum;
               else if (strcmp(name,"StressTensorxx")==0) return StressTensorxxEnum;
-              else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum;
-              else if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum;
          else stage=5;
    }
    if(stage==5){
-              if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum;
+              if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum;
+              else if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum;
+              else if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum;
               else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum;
               else if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum;

issm/trunk-jpl/src/m/enum/EnumDefinitions.py

@@ -330 +330 @@
 def GiaSolutionEnum(): return StringToEnum("GiaSolution")[0]
 def GiaAnalysisEnum(): return StringToEnum("GiaAnalysis")[0]
+def MeshdeformationSolutionEnum(): return StringToEnum("MeshdeformationSolution")[0]
+def MeshdeformationAnalysisEnum(): return StringToEnum("MeshdeformationAnalysis")[0]
 def ApproximationEnum(): return StringToEnum("Approximation")[0]
 def NoneApproximationEnum(): return StringToEnum("NoneApproximation")[0]