Index: /issm/trunk/src/c/objects/Icefront.cpp =================================================================== --- /issm/trunk/src/c/objects/Icefront.cpp (revision 3429) +++ /issm/trunk/src/c/objects/Icefront.cpp (revision 3430) @@ -16,5 +16,5 @@ #include "../include/typedefs.h" #include "../include/macros.h" -#include "./objects.h" +#include "./Icefront.h" @@ -25,6 +25,15 @@ } /*}}}*/ -/*FUNCTION Icefront creation {{{1*/ +/*FUNCTION Icefront constructor {{{1*/ Icefront::Icefront(char icefront_type[ICEFRONTSTRING],int icefront_fill,int icefront_sid, int icefront_mparid, int icefront_eid, int icefront_element_type, + int icefront_node_ids[MAX_ICEFRONT_GRIDS],double icefront_h[MAX_ICEFRONT_GRIDS],double icefront_b[MAX_ICEFRONT_GRIDS]){ + + this->Init(icefront_type,icefront_fill,icefront_sid, icefront_mparid, icefront_eid, icefront_element_type, + icefront_node_ids,icefront_h,icefront_b); + +} +/*}}}*/ +/*FUNCTION Init: used by constructor {{{1*/ +void Icefront::Init(char icefront_type[ICEFRONTSTRING],int icefront_fill,int icefront_sid, int icefront_mparid, int icefront_eid, int icefront_element_type, int icefront_node_ids[MAX_ICEFRONT_GRIDS],double icefront_h[MAX_ICEFRONT_GRIDS],double icefront_b[MAX_ICEFRONT_GRIDS]){ @@ -52,7 +61,79 @@ element=NULL; element_offset=UNDEF; - - return; -} +} +/*}}}*/ +/*FUNCTION Icefront constructor from iomodel {{{1*/ +Icefront::Icefront(int i, IoModel* iomodel){ + + int segment_width; + int element_type; + int i1,i2,i3,i4; + int element; + + /*icefront intermediary data: */ + char icefront_type[ICEFRONTSTRING]; + int icefront_fill; + int icefront_element_type; + int icefront_sid; + int icefront_eid; + int icefront_mparid; + int icefront_node_ids[MAX_ICEFRONT_GRIDS]; + double icefront_h[MAX_ICEFRONT_GRIDS]; + double icefront_b[MAX_ICEFRONT_GRIDS]; + + /*First, retrieve element index and element type: */ + if (strcmp(iomodel->meshtype,"2d")==0){ + segment_width=4; + element_type=TriaEnum(); + } + else{ + segment_width=6; + element_type=PentaEnum(); + } + element=(int)(*(iomodel->pressureload+segment_width*i+segment_width-2)-1); //element is in the penultimate column (grid1 grid2 ... elem fill) + + icefront_mparid=iomodel->numberofelements+1; //matlab indexing + icefront_mparid=iomodel->numberofelements+1; //matlab indexing + icefront_sid=i+1; //matlab indexing + icefront_fill=(int)*(iomodel->pressureload+segment_width*i+segment_width-1); + icefront_eid=(int) *(iomodel->pressureload+segment_width*i+segment_width-2); //matlab indexing + icefront_element_type=element_type; + + i1=(int)*(iomodel->pressureload+segment_width*i+0); + i2=(int)*(iomodel->pressureload+segment_width*i+1); + + icefront_node_ids[0]=i1; + icefront_node_ids[1]=i2; + + icefront_h[0]=iomodel->thickness[i1-1]; + icefront_h[1]=iomodel->thickness[i2-1]; + + icefront_b[0]=iomodel->bed[i1-1]; + icefront_b[1]=iomodel->bed[i2-1]; + + if ((int)*(iomodel->elements_type+2*element+0)==MacAyealFormulationEnum()){ //this is a collapsed 3d element, icefront will be 2d + strcpy(icefront_type,"segment"); + } + else if ((int)*(iomodel->elements_type+2*element+0)==PattynFormulationEnum()){ //this is a real 3d element, icefront will be 3d. + strcpy(icefront_type,"quad"); + i3=(int)*(iomodel->pressureload+segment_width*i+2); + i4=(int)*(iomodel->pressureload+segment_width*i+3); + icefront_node_ids[2]=i3; + icefront_node_ids[3]=i4; + + icefront_h[2]=iomodel->thickness[i3-1]; + icefront_h[3]=iomodel->thickness[i4-1]; + + icefront_b[2]=iomodel->bed[i3-1]; + icefront_b[3]=iomodel->bed[i4-1]; + } + else{ + ISSMERROR(exprintf(" element type %i not supported yet",(int)*(iomodel->elements_type+2*element+0))); + } + + this->Init(icefront_type,icefront_fill,icefront_sid,icefront_mparid,icefront_eid,icefront_element_type,icefront_node_ids,icefront_h,icefront_b); +} + + /*}}}*/ /*FUNCTION Icefront destructor {{{1*/ Index: /issm/trunk/src/c/objects/Icefront.h =================================================================== --- /issm/trunk/src/c/objects/Icefront.h (revision 3429) +++ /issm/trunk/src/c/objects/Icefront.h (revision 3430) @@ -16,4 +16,5 @@ #include "./Element.h" #include "./Node.h" +#include "../ModelProcessorx/IoModel.h" #define MAX_ICEFRONT_GRIDS 4 //max number of grids for a certain load @@ -23,7 +24,7 @@ private: - char type[ICEFRONTSTRING]; + char type[ICEFRONTSTRING]; int fill; - int sid; + int sid; /*material: */ @@ -49,22 +50,28 @@ public: + /*constructors: {{{1*/ Icefront(); Icefront(char type[ICEFRONTSTRING],int fill,int sid, int mparid, int eid, int element_type, int node_ids[MAX_ICEFRONT_GRIDS],double h[MAX_ICEFRONT_GRIDS],double b[MAX_ICEFRONT_GRIDS]); + void Init(char type[ICEFRONTSTRING],int fill,int sid, int mparid, int eid, int element_type, int node_ids[MAX_ICEFRONT_GRIDS],double h[MAX_ICEFRONT_GRIDS],double b[MAX_ICEFRONT_GRIDS]); + Icefront(int i, IoModel* iomodel); ~Icefront(); - + /*}}}*/ + /*object management: {{{1*/ + void Configure(void* elements,void* nodes,void* materials); + Object* copy(); + void DeepEcho(); + void Demarshall(char** pmarshalled_dataset); void Echo(); - void DeepEcho(); + int Enum(); + int GetId(); + char* GetName(); void Marshall(char** pmarshalled_dataset); int MarshallSize(); - char* GetName(); - void Demarshall(char** pmarshalled_dataset); - int Enum(); - int GetId(); int MyRank(); + /*}}}*/ + /*numerics: {{{1*/ void DistributeNumDofs(int* numdofspernode,int analysis_type,int sub_analysis_type); - void Configure(void* elements,void* nodes,void* materials); void CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type); void UpdateFromInputs(void* inputs); - void CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type); void CreatePVectorDiagnosticHoriz( Vec pg, void* inputs, int analysis_type,int sub_analysis_type); @@ -80,5 +87,5 @@ void PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type); void PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type); - Object* copy(); + /*}}}*/ }; Index: /issm/trunk/src/c/objects/Node.cpp =================================================================== --- /issm/trunk/src/c/objects/Node.cpp (revision 3429) +++ /issm/trunk/src/c/objects/Node.cpp (revision 3430) @@ -52,5 +52,5 @@ } /*}}}*/ -/*FUNCTION Node constructor from iomodel continuous Galerkin{{{2*/ +/*FUNCTION Node constructor from iomodel{{{2*/ Node::Node(int i, IoModel* iomodel){ //i is the node index @@ -98,5 +98,4 @@ /*set single point constraints.: */ - /*FROM DIAGNOSTICHORIZ*/ if (strcmp(iomodel->meshtype,"3d")==0){ /*We have a 3d mesh, we may have collapsed elements, hence dead grids. Freeze them out: */ @@ -112,11 +111,10 @@ } } - /*FROM DIAGNOSTICSTOKES*/ /*On a 3d mesh, in stokes formualtions, only stokes grids are free, the others are frozen: */ if (iomodel->borderstokes[i]){ //freeze everything except pressure - node->FreezeDof(1); - node->FreezeDof(2); - node->FreezeDof(3); + this->FreezeDof(1); + this->FreezeDof(2); + this->FreezeDof(3); } else if (iomodel->gridonstokes[i]==0){ @@ -125,57 +123,4 @@ } } - /*FROM DIAGNOSTICHUTTER*/ - if (!iomodel->gridonhutter[i]){ - for(k=1;k<=numdofs;k++){ - node->FreezeDof(k); - } - } -} -/*}}}*/ -/*FUNCTION Node constructor from iomodel discontinuous Galerkin{{{2*/ -Node::Node(int i,int j,IoModel* iomodel){ - /* i -> index of the vertex in C indexing - * j -> index of the node in C indexing*/ - - int numdofs; - int partitionborder; - int vertex_id; - int upper_node_id; - - /*id: */ - this->id=j+1; //matlab indexing - - /*indexing:*/ - DistributeNumDofs(&numdofs,iomodel->analysis_type,iomodel->sub_analysis_type); //number of dofs per node - if(iomodel->my_bordervertices[i])partitionborder=1; else partitionborder=0;//is this node on a partition border? - - this->indexing.Init(numdofs,partitionborder); - - /*properties (come from vertex number i): */ - this->properties.Init( - (int)iomodel->gridonbed[i], - (int)iomodel->gridonsurface[i], - (int)iomodel->gridoniceshelf[i], - (int)iomodel->gridonicesheet[i]); - - /*hooks: */ - vertex_id=i+1; //matlab indexing - - if (strcmp(iomodel->meshtype,"3d")==0){ - if (isnan(iomodel->uppernodes[i])){ - upper_node_id=this->id; //nodes on surface do not have upper nodes, only themselves. - } - else{ - upper_node_id=(int)iomodel->uppernodes[i]; - } - } - else{ - /*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/ - upper_node_id=this->id; - } - - this->hvertex.Init(&vertex_id,1); - this->hupper_node.Init(&upper_node_id,1); - } /*}}}*/ Index: /issm/trunk/src/c/objects/Node.h =================================================================== --- /issm/trunk/src/c/objects/Node.h (revision 3429) +++ /issm/trunk/src/c/objects/Node.h (revision 3430) @@ -44,5 +44,4 @@ Node(int id,DofIndexing* indexing, NodeProperties* properties, Hook* vertex, Hook* upper_node); Node(int i, IoModel* iomodel); - Node(int i,int j,IoModel* iomodel); ~Node(); /*}}}*/ Index: /issm/trunk/src/c/objects/Penpair.cpp =================================================================== --- /issm/trunk/src/c/objects/Penpair.cpp (revision 3429) +++ /issm/trunk/src/c/objects/Penpair.cpp (revision 3430) @@ -216,5 +216,6 @@ /*FUNCTION Penpair::PenaltyCreateKMatrix {{{1*/ void Penpair::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){ - + + /*If you code this piece, don't forget that a penalty will be inactive if it is dealing with clone nodes*/ /*No loads applied, do nothing: */ return; Index: /issm/trunk/src/c/objects/Riftfront.cpp =================================================================== --- /issm/trunk/src/c/objects/Riftfront.cpp (revision 3429) +++ /issm/trunk/src/c/objects/Riftfront.cpp (revision 3430) @@ -20,5 +20,5 @@ /*Object constructors and destructor*/ -/*FUNCTION Riftfront::constructor {{{1*/ +/*FUNCTION Riftfront::default constructor {{{1*/ Riftfront::Riftfront(){ /*in case :*/ @@ -27,6 +27,99 @@ } /*}}}1*/ -/*FUNCTION Riftfront::creation {{{1*/ +/*FUNCTION Riftfront::constructor {{{1*/ Riftfront::Riftfront(char riftfront_type[RIFTFRONTSTRING],int riftfront_id, int riftfront_node_ids[MAX_RIFTFRONT_GRIDS], int riftfront_mparid, double riftfront_h[MAX_RIFTFRONT_GRIDS],double riftfront_b[MAX_RIFTFRONT_GRIDS],double riftfront_s[MAX_RIFTFRONT_GRIDS],double riftfront_normal[2],double riftfront_length,int riftfront_fill,double riftfront_friction, double riftfront_fraction,double riftfront_fractionincrement, double riftfront_penalty_offset, int riftfront_penalty_lock, bool riftfront_active,bool riftfront_frozen, int riftfront_counter,bool riftfront_prestable,bool riftfront_shelf){ + + this->Init(riftfront_type,riftfront_id, riftfront_node_ids[MAX_RIFTFRONT_GRIDS], riftfront_mparid, riftfront_h[MAX_RIFTFRONT_GRIDS],riftfront_b[MAX_RIFTFRONT_GRIDS],riftfront_s[MAX_RIFTFRONT_GRIDS],riftfront_normal[2],riftfront_length,riftfront_fill,riftfront_friction, riftfront_fraction,riftfront_fractionincrement, riftfront_penalty_offset, riftfront_penalty_lock, griftfront_active,griftfront_frozen, riftfront_counter,riftfront_prestable,riftfront_shelf); + +} +/*}}}1*/ +/*FUNCTION Riftfront::constructor from iomodel{{{1*/ +Riftfront::Riftfront(int i, IoModel* iomodel){ + + /*rifts: */ + char riftfront_type[RIFTFRONTSTRING]; + int riftfront_id; + int riftfront_node_ids[2]; + int riftfront_mparid; + double riftfront_h[2]; + double riftfront_b[2]; + double riftfront_s[2]; + double riftfront_normal[2]; + double riftfront_length; + int riftfront_fill; + double riftfront_friction; + double riftfront_fraction; + double riftfront_fractionincrement; + bool riftfront_shelf; + double riftfront_penalty_offset; + int riftfront_penalty_lock; + bool riftfront_active; + bool riftfront_frozen; + int riftfront_counter; + bool riftfront_prestable; + int el1,el2; + int grid1,grid2; + double normal[2]; + double length; + int fill; + double friction; + double fraction; + double fractionincrement; + + /*Ok, retrieve all the data needed to add a penalty between the two grids: */ + el1=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+2); + el2=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+3); + + grid1=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+0); + grid2=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+1); + + normal[0]=*(iomodel->riftinfo+RIFTINFOSIZE*i+4); + normal[1]=*(iomodel->riftinfo+RIFTINFOSIZE*i+5); + length=*(iomodel->riftinfo+RIFTINFOSIZE*i+6); + + fill = (int)*(iomodel->riftinfo+RIFTINFOSIZE*i+7); + friction=*(iomodel->riftinfo+RIFTINFOSIZE*i+8); + fraction=*(iomodel->riftinfo+RIFTINFOSIZE*i+9); + fractionincrement=*(iomodel->riftinfo+RIFTINFOSIZE*i+10); + + strcpy(riftfront_type,"2d"); + riftfront_id=i+1; //matlab indexing + riftfront_node_ids[0]=grid1; + riftfront_node_ids[1]=grid2; + riftfront_mparid=iomodel->numberofelements+1; //matlab indexing + + riftfront_h[0]=iomodel->thickness[grid1-1]; + riftfront_h[1]=iomodel->thickness[grid2-1]; + + riftfront_b[0]=iomodel->bed[grid1-1]; + riftfront_b[1]=iomodel->bed[grid2-1]; + + riftfront_s[0]=iomodel->surface[grid1-1]; + riftfront_s[1]=iomodel->surface[grid2-1]; + + riftfront_normal[0]=normal[0]; + riftfront_normal[1]=normal[1]; + riftfront_length=length; + + riftfront_fill=fill; + riftfront_friction=friction; + riftfront_fraction=fraction; + riftfront_fractionincrement=fractionincrement; + riftfront_shelf=(bool)iomodel->gridoniceshelf[grid1-1]; + + riftfront_penalty_offset=iomodel->penalty_offset; + riftfront_penalty_lock=iomodel->penalty_lock; + + riftfront_active=0; + riftfront_frozen=0; + riftfront_counter=0; + riftfront_prestable=0; + + this->Init(riftfront_type,riftfront_id, riftfront_node_ids[MAX_RIFTFRONT_GRIDS], riftfront_mparid, riftfront_h[MAX_RIFTFRONT_GRIDS],riftfront_b[MAX_RIFTFRONT_GRIDS],riftfront_s[MAX_RIFTFRONT_GRIDS],riftfront_normal[2],riftfront_length,riftfront_fill,riftfront_friction, riftfront_fraction,riftfront_fractionincrement, riftfront_penalty_offset, riftfront_penalty_lock, griftfront_active,griftfront_frozen, riftfront_counter,riftfront_prestable,riftfront_shelf); + +} +/*}}}1*/ +/*FUNCTION Riftfront::Init, used by constructor {{{1*/ +void Riftfront::Init(char riftfront_type[RIFTFRONTSTRING],int riftfront_id, int riftfront_node_ids[MAX_RIFTFRONT_GRIDS], int riftfront_mparid, double riftfront_h[MAX_RIFTFRONT_GRIDS],double riftfront_b[MAX_RIFTFRONT_GRIDS],double riftfront_s[MAX_RIFTFRONT_GRIDS],double riftfront_normal[2],double riftfront_length,int riftfront_fill,double riftfront_friction, double riftfront_fraction,double riftfront_fractionincrement, double riftfront_penalty_offset, int riftfront_penalty_lock, bool riftfront_active,bool riftfront_frozen, int riftfront_counter,bool riftfront_prestable,bool riftfront_shelf){ int i; Index: /issm/trunk/src/c/objects/Riftfront.h =================================================================== --- /issm/trunk/src/c/objects/Riftfront.h (revision 3429) +++ /issm/trunk/src/c/objects/Riftfront.h (revision 3430) @@ -58,8 +58,11 @@ public: + /*constructors,destructors: {{{1*/ Riftfront(); + void Init(char type[RIFTFRONTSTRING],int id, int node_ids[MAX_RIFTFRONT_GRIDS], int mparid, double h[MAX_RIFTFRONT_GRIDS],double b[MAX_RIFTFRONT_GRIDS],double s[MAX_RIFTFRONT_GRIDS],double normal[2],double length,int fill,double friction, double fraction, double fractionincrement, double penalty_offset, int penalty_lock,bool active,bool frozen, int counter,bool prestable,bool shelf); Riftfront(char type[RIFTFRONTSTRING],int id, int node_ids[MAX_RIFTFRONT_GRIDS], int mparid, double h[MAX_RIFTFRONT_GRIDS],double b[MAX_RIFTFRONT_GRIDS],double s[MAX_RIFTFRONT_GRIDS],double normal[2],double length,int fill,double friction, double fraction, double fractionincrement, double penalty_offset, int penalty_lock,bool active,bool frozen, int counter,bool prestable,bool shelf); ~Riftfront(); - + /*}}}*/ + /*object management: {{{1*/ void Echo(); void DeepEcho(); @@ -72,4 +75,7 @@ int MyRank(); void Configure(void* elements,void* nodes,void* materials); + Object* copy(); + /*}}}*/ + /*numerics: {{{1*/ void UpdateFromInputs(void* inputs); void GetDofList(int* doflist,int* pnumberofdofs); @@ -89,5 +95,5 @@ int IsMaterialStable(void* inputs, int analysis_type); void OutputProperties(Vec riftproperties); - Object* copy(); + /*}}}*/ };