Index: /issm/trunk-jpl/m4/issm_options.m4 =================================================================== --- /issm/trunk-jpl/m4/issm_options.m4 (revision 17493) +++ /issm/trunk-jpl/m4/issm_options.m4 (revision 17494) @@ -1711,4 +1711,18 @@ AC_MSG_RESULT($HAVE_DAMAGEEVOLUTION) dnl }}} + dnl with-DepthAverage{{{ + AC_ARG_WITH([DepthAverage], + AS_HELP_STRING([--with-DepthAverage = YES], [compile with DepthAverage capabilities (default is yes)]), + [DEPTHAVERAGE=$withval],[DEPTHAVERAGE=yes]) + AC_MSG_CHECKING(for DepthAverage capability compilation) + + HAVE_DEPTHAVERAGE=no + if test "x$DEPTHAVERAGE" = "xyes"; then + HAVE_DEPTHAVERAGE=yes + AC_DEFINE([_HAVE_DEPTHAVERAGE_],[1],[with DepthAverage capability]) + fi + AM_CONDITIONAL([DEPTHAVERAGE], [test x$HAVE_DEPTHAVERAGE = xyes]) + AC_MSG_RESULT($HAVE_DEPTHAVERAGE) + dnl }}} dnl with-stressbalance{{{ AC_ARG_WITH([stressbalance], Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 17493) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 17494) @@ -368,4 +368,5 @@ ./cores/steadystate_core.cpp\ ./cores/masstransport_core.cpp\ + ./cores/depthaverage_core.cpp\ ./cores/extrudefrombase_core.cpp\ ./cores/extrudefromtop_core.cpp\ @@ -574,4 +575,7 @@ issm_sources += ./analyses/ExtrudeFromTopAnalysis.cpp endif +if DEPTHAVERAGE +issm_sources += ./analyses/DepthAverageAnalysis.cpp +endif if THERMAL issm_sources += ./analyses/ThermalAnalysis.cpp Index: /issm/trunk-jpl/src/c/analyses/DepthAverageAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/DepthAverageAnalysis.cpp (revision 17494) +++ /issm/trunk-jpl/src/c/analyses/DepthAverageAnalysis.cpp (revision 17494) @@ -0,0 +1,133 @@ +#include "./DepthAverageAnalysis.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" + +/*Model processing*/ +int DepthAverageAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ + return 1; +}/*}}}*/ +void DepthAverageAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ +}/*}}}*/ +void DepthAverageAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + + int counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum); + counter++; + } + } + + if(iomodel->meshtype==Mesh2DverticalEnum){ + iomodel->FetchDataToInput(elements,MeshVertexonbedEnum); + } +}/*}}}*/ +void DepthAverageAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ + + ::CreateNodes(nodes,iomodel,DepthAverageAnalysisEnum,P1Enum); + +}/*}}}*/ +void DepthAverageAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ +}/*}}}*/ +void DepthAverageAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ +}/*}}}*/ + +/*Finite Element Analysis*/ +void DepthAverageAnalysis::Core(FemModel* femmodel){/*{{{*/ + _error_("not implemented"); +}/*}}}*/ +ElementVector* DepthAverageAnalysis::CreateDVector(Element* element){/*{{{*/ + /*Default, return NULL*/ + return NULL; +}/*}}}*/ +ElementMatrix* DepthAverageAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/ +_error_("Not implemented"); +}/*}}}*/ +ElementMatrix* DepthAverageAnalysis::CreateKMatrix(Element* element){/*{{{*/ + + /*Intermediaries */ + IssmDouble Jdet,D; + IssmDouble *xyz_list = NULL; + + /*Get dimension*/ + int dim,meshtype; + element->FindParam(&meshtype,MeshTypeEnum); + switch(meshtype){ + case Mesh2DverticalEnum: dim = 2; break; + case Mesh3DEnum: dim = 3; break; + default: _error_("mesh "<GetNumberOfNodes(); + + /*Initialize Element vector and other vectors*/ + ElementMatrix* Ke = element->NewElementMatrix(); + IssmDouble* B = xNew(numnodes); + IssmDouble* Bprime = xNew(numnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(2); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctions(Bprime,gauss); + GetB(B,element,dim,xyz_list,gauss); + D=gauss->weight*Jdet; + + TripleMultiply(B,1,numnodes,1, + &D,1,1,0, + Bprime,1,numnodes,0, + &Ke->values[0],1); + } + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(B); + xDelete(Bprime); + delete gauss; + return Ke; +}/*}}}*/ +ElementVector* DepthAverageAnalysis::CreatePVector(Element* element){/*{{{*/ + return NULL; +}/*}}}*/ +void DepthAverageAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ + /* Compute B matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz]; + where hi is the interpolation function for node i.*/ + + /*Fetch number of nodes for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Get nodal functions derivatives*/ + IssmDouble* dbasis=xNew(dim*numnodes); + element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + + /*Build B: */ + for(int i=0;i(dbasis); +} +/*}}}*/ +void DepthAverageAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ +void DepthAverageAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ + + int inputenum; + element->FindParam(&inputenum,InputToDepthaverageEnum); + element->InputUpdateFromSolutionOneDof(solution,inputenum); +}/*}}}*/ +void DepthAverageAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ + /*Default, do nothing*/ + return; +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/DepthAverageAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/DepthAverageAnalysis.h (revision 17494) +++ /issm/trunk-jpl/src/c/analyses/DepthAverageAnalysis.h (revision 17494) @@ -0,0 +1,33 @@ +/*! \file DepthAverageAnalysis.h + * \brief: header file for generic external result object + */ + +#ifndef _DepthAverageAnalysis_ +#define _DepthAverageAnalysis_ + +/*Headers*/ +#include "./Analysis.h" + +class DepthAverageAnalysis: public Analysis{ + + public: + /*Model processing*/ + int DofsPerNode(int** doflist,int meshtype,int approximation); + void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void CreateNodes(Nodes* nodes,IoModel* iomodel); + void CreateConstraints(Constraints* constraints,IoModel* iomodel); + void CreateLoads(Loads* loads, IoModel* iomodel); + + /*Finite element Analysis*/ + void Core(FemModel* femmodel); + ElementVector* CreateDVector(Element* element); + ElementMatrix* CreateJacobianMatrix(Element* element); + ElementMatrix* CreateKMatrix(Element* element); + ElementVector* CreatePVector(Element* element); + void GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); + void GetSolutionFromInputs(Vector* solution,Element* element); + void InputUpdateFromSolution(IssmDouble* solution,Element* element); + void UpdateConstraints(FemModel* femmodel); +}; +#endif Index: /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp (revision 17493) +++ /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp (revision 17494) @@ -77,4 +77,7 @@ case ExtrudeFromTopAnalysisEnum : return new ExtrudeFromTopAnalysis(); #endif + #ifdef _HAVE_DEPTHAVERAGE_ + case DepthAverageAnalysisEnum : return new DepthAverageAnalysis(); + #endif #ifdef _HAVE_SMOOTHEDSURFACESLOPEX_ case SmoothedSurfaceSlopeXAnalysisEnum : return new SmoothedSurfaceSlopeXAnalysis(); Index: /issm/trunk-jpl/src/c/analyses/analyses.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/analyses.h (revision 17493) +++ /issm/trunk-jpl/src/c/analyses/analyses.h (revision 17494) @@ -14,4 +14,5 @@ #include "./BalancevelocityAnalysis.h" #include "./DamageEvolutionAnalysis.h" +#include "./DepthAverageAnalysis.h" #include "./EnthalpyAnalysis.h" #include "./ExtrudeFromBaseAnalysis.h" Index: /issm/trunk-jpl/src/c/cores/depthaverage_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/depthaverage_core.cpp (revision 17494) +++ /issm/trunk-jpl/src/c/cores/depthaverage_core.cpp (revision 17494) @@ -0,0 +1,20 @@ +/*!\file: depthaverage_core.cpp + * \brief: core of the extrusion solution + */ + +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../solutionsequences/solutionsequences.h" +#include "../modules/modules.h" + +void depthaverage_core(FemModel* femmodel){ + + if(VerboseSolution()) _printf0_("extruding solution from base...\n"); + + /*Call on core computations: */ + femmodel->SetCurrentConfiguration(DepthAverageAnalysisEnum); + solutionsequence_linear(femmodel); + +} Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 17493) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 17494) @@ -346,4 +346,5 @@ ExtrudeFromBaseAnalysisEnum, ExtrudeFromTopAnalysisEnum, + DepthAverageAnalysisEnum, SteadystateSolutionEnum, SurfaceSlopeSolutionEnum, @@ -406,4 +407,5 @@ InputToExtrudeEnum, InputToL2ProjectEnum, + InputToDepthaverageEnum, IntParamEnum, IntVecParamEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 17493) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 17494) @@ -347,4 +347,5 @@ case ExtrudeFromBaseAnalysisEnum : return "ExtrudeFromBaseAnalysis"; case ExtrudeFromTopAnalysisEnum : return "ExtrudeFromTopAnalysis"; + case DepthAverageAnalysisEnum : return "DepthAverageAnalysis"; case SteadystateSolutionEnum : return "SteadystateSolution"; case SurfaceSlopeSolutionEnum : return "SurfaceSlopeSolution"; @@ -401,4 +402,5 @@ case InputToExtrudeEnum : return "InputToExtrude"; case InputToL2ProjectEnum : return "InputToL2Project"; + case InputToDepthaverageEnum : return "InputToDepthaverage"; case IntParamEnum : return "IntParam"; case IntVecParamEnum : return "IntVecParam"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 17493) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 17494) @@ -353,4 +353,5 @@ else if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum; else if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum; + else if (strcmp(name,"DepthAverageAnalysis")==0) return DepthAverageAnalysisEnum; else if (strcmp(name,"SteadystateSolution")==0) return SteadystateSolutionEnum; else if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum; @@ -382,9 +383,9 @@ else if (strcmp(name,"Loads")==0) return LoadsEnum; else if (strcmp(name,"Materials")==0) return MaterialsEnum; - else if (strcmp(name,"Nodes")==0) return NodesEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"Contours")==0) return ContoursEnum; + if (strcmp(name,"Nodes")==0) return NodesEnum; + else if (strcmp(name,"Contours")==0) return ContoursEnum; else if (strcmp(name,"Parameters")==0) return ParametersEnum; else if (strcmp(name,"Vertices")==0) return VerticesEnum; @@ -410,4 +411,5 @@ else if (strcmp(name,"InputToExtrude")==0) return InputToExtrudeEnum; else if (strcmp(name,"InputToL2Project")==0) return InputToL2ProjectEnum; + else if (strcmp(name,"InputToDepthaverage")==0) return InputToDepthaverageEnum; else if (strcmp(name,"IntParam")==0) return IntParamEnum; else if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum; @@ -504,10 +506,10 @@ else if (strcmp(name,"Vy")==0) return VyEnum; else if (strcmp(name,"VyPicard")==0) return VyPicardEnum; - else if (strcmp(name,"Vz")==0) return VzEnum; - else if (strcmp(name,"VzSSA")==0) return VzSSAEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"VzHO")==0) return VzHOEnum; + if (strcmp(name,"Vz")==0) return VzEnum; + else if (strcmp(name,"VzSSA")==0) return VzSSAEnum; + else if (strcmp(name,"VzHO")==0) return VzHOEnum; else if (strcmp(name,"VzPicard")==0) return VzPicardEnum; else if (strcmp(name,"VzFS")==0) return VzFSEnum; @@ -627,10 +629,10 @@ else if (strcmp(name,"OutputFilePointer")==0) return OutputFilePointerEnum; else if (strcmp(name,"OutputFileName")==0) return OutputFileNameEnum; - else if (strcmp(name,"LockFileName")==0) return LockFileNameEnum; - else if (strcmp(name,"ToolkitsOptionsAnalyses")==0) return ToolkitsOptionsAnalysesEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"ToolkitsOptionsStrings")==0) return ToolkitsOptionsStringsEnum; + if (strcmp(name,"LockFileName")==0) return LockFileNameEnum; + else if (strcmp(name,"ToolkitsOptionsAnalyses")==0) return ToolkitsOptionsAnalysesEnum; + else if (strcmp(name,"ToolkitsOptionsStrings")==0) return ToolkitsOptionsStringsEnum; else if (strcmp(name,"QmuErrName")==0) return QmuErrNameEnum; else if (strcmp(name,"QmuInName")==0) return QmuInNameEnum; Index: /issm/trunk-jpl/src/c/shared/Exceptions/Exceptions.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Exceptions/Exceptions.cpp (revision 17493) +++ /issm/trunk-jpl/src/c/shared/Exceptions/Exceptions.cpp (revision 17494) @@ -46,5 +46,5 @@ file_line= what_line; /*When error messages are not shown properly, uncomment the following line*/ - //this->Report(); + this->Report(); }/*}}}*/