Index: sm/trunk-jpl/src/c/classes/Loads/Icefront.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Icefront.cpp (revision 15648) +++ (revision ) @@ -1,932 +1,0 @@ -/*!\file Icefront.c - * \brief: implementation of the Icefront object - */ - -/*Headers:*/ -/*{{{*/ -#ifdef HAVE_CONFIG_H - #include -#else -#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" -#endif - -#include "../classes.h" -#include "shared/shared.h" -/*}}}*/ - -/*Load macros*/ -#define NUMVERTICESSEG 2 -#define NUMVERTICESQUA 4 - -/*Icefront constructors and destructor*/ -/*FUNCTION Icefront::Icefront() {{{*/ -Icefront::Icefront(){ - - this->inputs=NULL; - this->parameters=NULL; - - this->hnodes=NULL; - this->nodes= NULL; - this->hvertices=NULL; - this->vertices= NULL; - this->helement=NULL; - this->element= NULL; - this->hmatpar=NULL; - this->matpar= NULL; -} -/*}}}*/ -/*FUNCTION Icefront::Icefront(int id, int i, IoModel* iomodel,int analysis_type) {{{*/ -Icefront::Icefront(int icefront_id,int i, IoModel* iomodel,int in_icefront_type, int in_analysis_type){ - - int segment_width; - int element; - int numnodes; - int numvertices; - - /*icefront constructor data: */ - int icefront_eid; - int icefront_mparid; - int icefront_fill; - int *icefront_node_ids = NULL; - int *icefront_vertex_ids = NULL; - - /*First, retrieve element index and element type: */ - if(iomodel->dim==2){ - segment_width=4; - } - else{ - segment_width=6; - } - _assert_(iomodel->Data(DiagnosticIcefrontEnum)); - element=reCast(iomodel->Data(DiagnosticIcefrontEnum)[segment_width*i+segment_width-2]-1); - - /*Build ids for hook constructors: */ - icefront_eid = reCast(iomodel->Data(DiagnosticIcefrontEnum)[segment_width*i+segment_width-2]); //matlab indexing - icefront_mparid = iomodel->numberofelements+1; - - if (in_icefront_type==SSA2dIceFrontEnum || in_icefront_type==SSA3dIceFrontEnum){ - numnodes = 2; - numvertices = 2; - icefront_node_ids = xNew(numnodes); - icefront_vertex_ids = xNew(numvertices); - icefront_node_ids[0]=iomodel->nodecounter+reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0)); - icefront_node_ids[1]=iomodel->nodecounter+reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1)); - icefront_vertex_ids[0]=reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0)); - icefront_vertex_ids[1]=reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1)); - } - else if (in_icefront_type==HOIceFrontEnum || in_icefront_type==FSIceFrontEnum){ - numnodes = 4; - numvertices = 4; - icefront_node_ids = xNew(numnodes); - icefront_vertex_ids = xNew(numvertices); - icefront_node_ids[0]=iomodel->nodecounter+reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0)); - icefront_node_ids[1]=iomodel->nodecounter+reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1)); - icefront_node_ids[2]=iomodel->nodecounter+reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+2)); - icefront_node_ids[3]=iomodel->nodecounter+reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+3)); - icefront_vertex_ids[0]=reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0)); - icefront_vertex_ids[1]=reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1)); - icefront_vertex_ids[2]=reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+2)); - icefront_vertex_ids[3]=reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+3)); - } - else _error_("in_icefront_type " << EnumToStringx(in_icefront_type) << " not supported yet!"); - - /*Fill*/ - icefront_fill=reCast(iomodel->Data(DiagnosticIcefrontEnum)[segment_width*i+segment_width-1]); - - /*Ok, we have everything to build the object: */ - this->id=icefront_id; - this->analysis_type=in_analysis_type; - - /*Hooks: */ - this->hnodes=new Hook(icefront_node_ids,numnodes); - this->hvertices=new Hook(icefront_vertex_ids,numvertices); - this->helement=new Hook(&icefront_eid,1); - this->hmatpar=new Hook(&icefront_mparid,1); - - //intialize and add as many inputs per element as requested: - this->inputs=new Inputs(); - this->inputs->AddInput(new IntInput(FillEnum,icefront_fill)); - this->inputs->AddInput(new IntInput(IceFrontTypeEnum,in_icefront_type)); - - //parameters and hooked fields: we still can't point to them, they may not even exist. Configure will handle this. - this->parameters = NULL; - this->nodes = NULL; - this->vertices = NULL; - this->element = NULL; - this->matpar = NULL; - - /*Clean up*/ - xDelete(icefront_node_ids); - xDelete(icefront_vertex_ids); -} -/*}}}*/ -/*FUNCTION Icefront::Icefront(const char* element_type_in,Inputs* inputs_in,Matpar* matpar_in, int icefront_type, int in_analysis_type) {{{*/ -Icefront::Icefront(int element_id, const char* element_type_in,Inputs* inputs_in,Matpar* matpar_in,int in_icefront_type, int in_analysis_type){ - - int segment_width; - int element; - int numnodes; - int numvertices; - int numberofelements; - - /*icefront constructor data: */ - int icefront_eid; - int icefront_mparid; - int icefront_node_ids[NUMVERTICESQUA]; //initialize with largest size - int icefront_vertex_ids[NUMVERTICESQUA]; //initialize with largest size - -// /*find parameters: */ -// iomodel->Constant(&numberofelements,MeshNumberofelementsEnum); -// - /*First, retrieve element index and element type: */ - if(strcmp(element_type_in,"2d")==0){ - segment_width=4; - } - else{ - segment_width=6; - } -// element=element_number; - - /*Build ids for hook constructors: */ -// icefront_eid=reCast( *(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+segment_width-2)); //matlab indexing -// icefront_mparid=numberofelements+1; //matlab indexing -// - if (in_icefront_type==SSA2dIceFrontEnum || in_icefront_type==SSA3dIceFrontEnum){ -// icefront_node_ids[0]=iomodel->nodecounter+reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0)); -// icefront_node_ids[1]=iomodel->nodecounter+reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1)); -// icefront_vertex_ids[0]=reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0)); -// icefront_vertex_ids[1]=reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1)); - } - else if (in_icefront_type==HOIceFrontEnum || in_icefront_type==FSIceFrontEnum){ -// icefront_node_ids[0]=iomodel->nodecounter+reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0)); -// icefront_node_ids[1]=iomodel->nodecounter+reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1)); -// icefront_node_ids[2]=iomodel->nodecounter+reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+2)); -// icefront_node_ids[3]=iomodel->nodecounter+reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+3)); -// icefront_vertex_ids[0]=reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0)); -// icefront_vertex_ids[1]=reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1)); -// icefront_vertex_ids[2]=reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+2)); -// icefront_vertex_ids[3]=reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+3)); - } - else _error_("in_icefront_type " << EnumToStringx(in_icefront_type) << " not supported yet!"); - - if (in_icefront_type==HOIceFrontEnum || in_icefront_type==FSIceFrontEnum){ - numnodes=4; - numvertices=4; - } - else{ - numnodes=2; - numvertices=2; - } - - /*Ok, we have everything to build the object: */ - this->id=1; - this->analysis_type=in_analysis_type; - - /*Hooks: */ - this->hnodes=new Hook(icefront_node_ids,numnodes); - this->hvertices=new Hook(icefront_vertex_ids,numvertices); - this->helement=new Hook(&icefront_eid,1); - this->hmatpar=new Hook(&icefront_mparid,1); - - //intialize and add as many inputs per element as requested: - this->inputs=inputs_in; - this->inputs->AddInput(new IntInput(FillEnum,1)); //We always consider we have water, if above sea level, only air will be applied - this->inputs->AddInput(new IntInput(IceFrontTypeEnum,in_icefront_type)); - - //parameters and hooked fields: we still can't point to them, they may not even exist. Configure will handle this. - this->parameters = NULL; - this->nodes = NULL; - this->vertices = NULL; - this->element = NULL; - this->matpar = matpar_in; -} - -/*}}}*/ -/*FUNCTION Icefront::~Icefront() {{{*/ -Icefront::~Icefront(){ - delete inputs; - this->parameters=NULL; - delete hnodes; - delete hvertices; - delete helement; - delete hmatpar; -} -/*}}}*/ - -/*Object virtual functions definitions:*/ -/*FUNCTION Icefront::Echo {{{*/ -void Icefront::Echo(void){ - _printf_("Icefront:\n"); - _printf_(" id: " << id << "\n"); - _printf_(" analysis_type: " << EnumToStringx(analysis_type) << "\n"); - hnodes->Echo(); - hvertices->Echo(); - helement->Echo(); - hmatpar->Echo(); - _printf_(" parameters: " << parameters << "\n"); - if(parameters)parameters->Echo(); - _printf_(" inputs: " << inputs << "\n"); - if(inputs)inputs->Echo(); -} -/*}}}*/ -/*FUNCTION Icefront::DeepEcho{{{*/ -void Icefront::DeepEcho(void){ - - _printf_("Icefront:\n"); - _printf_(" id: " << id << "\n"); - _printf_(" analysis_type: " << EnumToStringx(analysis_type) << "\n"); - hnodes->DeepEcho(); - hvertices->DeepEcho(); - helement->DeepEcho(); - hmatpar->DeepEcho(); - _printf_(" parameters: " << parameters << "\n"); - if(parameters)parameters->DeepEcho(); - _printf_(" inputs: " << inputs << "\n"); - if(inputs)inputs->DeepEcho(); -} -/*}}}*/ -/*FUNCTION Icefront::Id {{{*/ -int Icefront::Id(void){ return id; } -/*}}}*/ -/*FUNCTION Icefront::ObjectEnum{{{*/ -int Icefront::ObjectEnum(void){ - - return IcefrontEnum; - -} -/*}}}*/ -/*FUNCTION Icefront::copy {{{*/ -Object* Icefront::copy() { - - Icefront* icefront=NULL; - - icefront=new Icefront(); - - /*copy fields: */ - icefront->id=this->id; - icefront->analysis_type=this->analysis_type; - if(this->inputs){ - icefront->inputs=(Inputs*)this->inputs->Copy(); - } - else{ - icefront->inputs=new Inputs(); - } - /*point parameters: */ - icefront->parameters=this->parameters; - - /*now deal with hooks and objects: */ - icefront->hnodes = (Hook*)this->hnodes->copy(); - icefront->hvertices = (Hook*)this->hvertices->copy(); - icefront->helement = (Hook*)this->helement->copy(); - icefront->hmatpar = (Hook*)this->hmatpar->copy(); - - /*corresponding fields*/ - icefront->nodes = (Node**)icefront->hnodes->deliverp(); - icefront->vertices = (Vertex**)icefront->hvertices->deliverp(); - icefront->element = (Element*)icefront->helement->delivers(); - icefront->matpar = (Matpar*)icefront->hmatpar->delivers(); - - return icefront; - -} -/*}}}*/ - -/*Load virtual functions definitions:*/ -/*FUNCTION Icefront::Configure {{{*/ -void Icefront::Configure(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){ - - /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective - * datasets, using internal ids and offsets hidden in hooks: */ - hnodes->configure((DataSet*)nodesin); - hvertices->configure((DataSet*)verticesin); - helement->configure((DataSet*)elementsin); - hmatpar->configure((DataSet*)materialsin); - - /*Initialize hooked fields*/ - this->nodes = (Node**)hnodes->deliverp(); - this->vertices = (Vertex**)hvertices->deliverp(); - this->element = (Element*)helement->delivers(); - this->matpar = (Matpar*)hmatpar->delivers(); - - /*point parameters to real dataset: */ - this->parameters=parametersin; -} -/*}}}*/ -/*FUNCTION Icefront::SetCurrentConfiguration {{{*/ -void Icefront::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){ -} -/*}}}*/ -/*FUNCTION Icefront::CreateKMatrix {{{*/ -void Icefront::CreateKMatrix(Matrix* Kff, Matrix* Kfs){ - - /*No stiffness loads applied, do nothing: */ - return; - -} -/*}}}*/ -/*FUNCTION Icefront::CreatePVector {{{*/ -void Icefront::CreatePVector(Vector* pf){ - - /*Checks in debugging mode*/ - _assert_(nodes); - _assert_(element); - _assert_(matpar); - - /*Retrieve parameters: */ - ElementVector* pe=NULL; - int analysis_type; - this->parameters->FindParam(&analysis_type,AnalysisTypeEnum); - - /*Just branch to the correct element icefront vector generator, according to the type of analysis we are carrying out: */ - switch(analysis_type){ - #ifdef _HAVE_DIAGNOSTIC_ - case DiagnosticHorizAnalysisEnum: - pe=CreatePVectorDiagnosticHoriz(); - break; - #endif - #ifdef _HAVE_CONTROL_ - case AdjointHorizAnalysisEnum: - pe=CreatePVectorAdjointHoriz(); - break; - #endif - default: - _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet"); - } - - /*Add to global Vector*/ - if(pe){ - pe->AddToGlobal(pf); - delete pe; - } -} -/*}}}*/ -/*FUNCTION Icefront::CreateJacobianMatrix{{{*/ -void Icefront::CreateJacobianMatrix(Matrix* Jff){ - this->CreateKMatrix(Jff,NULL); -} -/*}}}*/ -/*FUNCTION Icefront::GetNodesSidList{{{*/ -void Icefront::GetNodesSidList(int* sidlist){ - - int type; - inputs->GetInputValue(&type,IceFrontTypeEnum); - _assert_(sidlist); - _assert_(nodes); - - switch(type){ - case SSA2dIceFrontEnum: - case SSA3dIceFrontEnum: - for(int i=0;iSid(); - return; -#ifdef _HAVE_3D_ - case HOIceFrontEnum: - case FSIceFrontEnum: - for(int i=0;iSid(); - return; -#endif - default: - _error_("Icefront type " << EnumToStringx(type) << " not supported yet"); - } -} -/*}}}*/ -/*FUNCTION Icefront::GetNumberOfNodes{{{*/ -int Icefront::GetNumberOfNodes(void){ - - int type; - inputs->GetInputValue(&type,IceFrontTypeEnum); - - switch(type){ - case SSA2dIceFrontEnum: - return NUMVERTICESSEG; -#ifdef _HAVE_3D_ - case SSA3dIceFrontEnum: - return NUMVERTICESSEG; - case HOIceFrontEnum: - return NUMVERTICESQUA; - case FSIceFrontEnum: - return NUMVERTICESQUA; -#endif - default: - _error_("Icefront type " << EnumToStringx(type) << " not supported yet"); - } - -} -/*}}}*/ -/*FUNCTION Icefront::IsPenalty{{{*/ -bool Icefront::IsPenalty(void){ - return false; -} -/*}}}*/ -/*FUNCTION Icefront::PenaltyCreateKMatrix {{{*/ -void Icefront::PenaltyCreateKMatrix(Matrix* Kff, Matrix* Kfs, IssmDouble kmax){ - /*do nothing: */ - return; -} -/*}}}*/ -/*FUNCTION Icefront::PenaltyCreatePVector{{{*/ -void Icefront::PenaltyCreatePVector(Vector* pf,IssmDouble kmax){ - /*do nothing: */ - return; -} -/*}}}*/ -/*FUNCTION Icefront::PenaltyCreateJacobianMatrix{{{*/ -void Icefront::PenaltyCreateJacobianMatrix(Matrix* Jff,IssmDouble kmax){ - this->PenaltyCreateKMatrix(Jff,NULL,kmax); -} -/*}}}*/ -/*FUNCTION Icefront::SetwiseNodeConnectivity{{{*/ -void Icefront::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int set1_enum,int set2_enum){ - - /*Output */ - int d_nz = 0; - int o_nz = 0; - - /*Loop over all nodes*/ - for(int i=0;iGetNumberOfNodes();i++){ - - if(!flags[this->nodes[i]->Sid()]){ - - /*flag current node so that no other element processes it*/ - flags[this->nodes[i]->Sid()]=true; - - /*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/ - switch(set2_enum){ - case FsetEnum: - if(nodes[i]->indexing.fsize){ - if(this->nodes[i]->IsClone()) - o_nz += 1; - else - d_nz += 1; - } - break; - case GsetEnum: - if(nodes[i]->indexing.gsize){ - if(this->nodes[i]->IsClone()) - o_nz += 1; - else - d_nz += 1; - } - break; - case SsetEnum: - if(nodes[i]->indexing.ssize){ - if(this->nodes[i]->IsClone()) - o_nz += 1; - else - d_nz += 1; - } - break; - default: _error_("not supported"); - } - } - } - - /*Assign output pointers: */ - *pd_nz=d_nz; - *po_nz=o_nz; -} -/*}}}*/ -/*FUNCTION Icefront::InAnalysis{{{*/ -bool Icefront::InAnalysis(int in_analysis_type){ - if (in_analysis_type==this->analysis_type)return true; - else return false; -} -/*}}}*/ - -/*Update virtual functions definitions:*/ -/*FUNCTION Icefront::InputUpdateFromVector(IssmDouble* vector, int name, int type) {{{*/ -void Icefront::InputUpdateFromVector(IssmDouble* vector, int name, int type){ - /*Nothing updated yet*/ -} -/*}}}*/ -/*FUNCTION Icefront::InputUpdateFromVector(int* vector, int name, int type) {{{*/ -void Icefront::InputUpdateFromVector(int* vector, int name, int type){ - /*Nothing updated yet*/ -} -/*}}}*/ -/*FUNCTION Icefront::InputUpdateFromVector(bool* vector, int name, int type) {{{*/ -void Icefront::InputUpdateFromVector(bool* vector, int name, int type){ - /*Nothing updated yet*/ -} -/*}}}*/ -/*FUNCTION Icefront::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type) {{{*/ -void Icefront::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){ - /*Nothing updated yet*/ -} -/*}}}*/ -/*FUNCTION Icefront::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type) {{{*/ -void Icefront::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){ - /*Nothing updated yet*/ -} -/*}}}*/ -/*FUNCTION Icefront::InputUpdateFromVectorDakota(int* vector, int name, int type) {{{*/ -void Icefront::InputUpdateFromVectorDakota(int* vector, int name, int type){ - /*Nothing updated yet*/ -} -/*}}}*/ -/*FUNCTION Icefront::InputUpdateFromVectorDakota(bool* vector, int name, int type) {{{*/ -void Icefront::InputUpdateFromVectorDakota(bool* vector, int name, int type){ - /*Nothing updated yet*/ -} -/*}}}*/ -/*FUNCTION Icefront::InputUpdateFromConstant(IssmDouble constant, int name) {{{*/ -void Icefront::InputUpdateFromConstant(IssmDouble constant, int name){ - /*Nothing updated yet*/ -} -/*}}}*/ -/*FUNCTION Icefront::InputUpdateFromConstant(int constant, int name) {{{*/ -void Icefront::InputUpdateFromConstant(int constant, int name){ - /*Nothing updated yet*/ -} -/*}}}*/ -/*FUNCTION Icefront::InputUpdateFromConstant(bool constant, int name) {{{*/ -void Icefront::InputUpdateFromConstant(bool constant, int name){ - /*Nothing updated yet*/ -} -/*}}}*/ -/*FUNCTION Icefront::InputUpdateFromSolution{{{*/ -void Icefront::InputUpdateFromSolution(IssmDouble* solution){ - /*Nothing updated yet*/ -} -/*}}}*/ - -/*Icefront numerics: */ -#ifdef _HAVE_DIAGNOSTIC_ -/*FUNCTION Icefront::CreatePVectorDiagnosticHoriz {{{*/ -ElementVector* Icefront::CreatePVectorDiagnosticHoriz(void){ - - int type; - inputs->GetInputValue(&type,IceFrontTypeEnum); - - switch(type){ - case SSA2dIceFrontEnum: - return CreatePVectorDiagnosticSSA2d(); - #ifdef _HAVE_3D_ - case SSA3dIceFrontEnum: - return CreatePVectorDiagnosticSSA3d(); - case HOIceFrontEnum: - return CreatePVectorDiagnosticHO(); - case FSIceFrontEnum: - return CreatePVectorDiagnosticFS(); - #endif - default: - _error_("Icefront type " << EnumToStringx(type) << " not supported yet"); - } -} -/*}}}*/ -/*FUNCTION Icefront::CreatePVectorDiagnosticSSA2d{{{*/ -ElementVector* Icefront::CreatePVectorDiagnosticSSA2d(void){ - - /*Intermediary*/ - int ig,index1,index2,fill; - IssmDouble Jdet; - IssmDouble thickness,bed,pressure,ice_pressure,rho_water,rho_ice,gravity; - IssmDouble water_pressure,air_pressure,surface_under_water,base_under_water; - IssmDouble xyz_list[NUMVERTICESSEG][3]; - IssmDouble normal[2]; - GaussTria *gauss; - - /*return of element is on water*/ - Tria* tria=((Tria*)element); - if(tria->IsOnWater()) return NULL; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = tria->NumberofNodes(); - int numdof = numnodes*NDOF2; - - /*Initialize Element vector and vectors*/ - ElementVector* pe=new ElementVector(tria->nodes,numnodes,this->parameters,SSAApproximationEnum); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICESSEG); - Input* thickness_input=tria->inputs->GetInput(ThicknessEnum); _assert_(thickness_input); - Input* bed_input =tria->inputs->GetInput(BedEnum); _assert_(bed_input); - inputs->GetInputValue(&fill,FillEnum); - rho_water=matpar->GetRhoWater(); - rho_ice =matpar->GetRhoIce(); - gravity =matpar->GetG(); - GetSegmentNormal(&normal[0],xyz_list); - - /*Start looping on Gaussian points*/ - index1=tria->GetNodeIndex(nodes[0]); - index2=tria->GetNodeIndex(nodes[1]); - gauss=new GaussTria(index1,index2,3); - - for(ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - thickness_input->GetInputValue(&thickness,gauss); - bed_input->GetInputValue(&bed,gauss); - - switch(fill){ - case WaterEnum: - surface_under_water=min(0.,thickness+bed); // 0 if the top of the glacier is above water level - base_under_water=min(0.,bed); // 0 if the bottom of the glacier is above water level - water_pressure=1.0/2.0*gravity*rho_water*(pow(surface_under_water,2) - pow(base_under_water,2)); - break; - case AirEnum: - water_pressure=0; - break; - case IceEnum: - water_pressure=-1.0/2.0*gravity*rho_ice*pow(thickness,2); // we are facing a wall of ice. use water_pressure to cancel the lithostatic pressure. - break; - default: - _error_("fill type " << EnumToStringx(fill) << " not supported yet"); - } - ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2); - air_pressure=0; - pressure = ice_pressure + water_pressure + air_pressure; - - tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss); - tria->GetNodalFunctions(basis,gauss); - - for (int i=0;ivalues[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i]; - pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i]; - } - } - - /*Transform load vector*/ - TransformLoadVectorCoord(pe,tria->nodes,numnodes,XYEnum); - - /*Clean up and return*/ - xDelete(basis); - delete gauss; - return pe; -} -/*}}}*/ -#endif - -#ifdef _HAVE_CONTROL_ -/*FUNCTION Icefront::CreatePVectorAdjointHoriz {{{*/ -ElementVector* Icefront::CreatePVectorAdjointHoriz(void){ - - /*No load vector applied to the adjoint*/ - return NULL; -} -/*}}}*/ -#endif -#ifdef _HAVE_3D_ -/*FUNCTION Icefront::CreatePVectorDiagnosticSSA3d{{{*/ -ElementVector* Icefront::CreatePVectorDiagnosticSSA3d(void){ - - Icefront *icefront = NULL; - Penta *penta = NULL; - Tria *tria = NULL; - - /*Cast element onto Penta*/ - penta =(Penta*)this->element; - - /*Return if not on bed*/ - if(!penta->IsOnBed() || penta->IsOnWater()) return NULL; - - /*Spawn Tria and call SSA2d*/ - tria =(Tria*)penta->SpawnTria(0,1,2); - icefront=(Icefront*)this->copy(); - icefront->element=tria; - icefront->inputs->AddInput(new IntInput(IceFrontTypeEnum,SSA2dIceFrontEnum)); - ElementVector* pe=icefront->CreatePVectorDiagnosticSSA2d(); - - /*clean-up and return*/ - delete tria->material; - delete tria; - delete icefront; - return pe; -} -/*}}}*/ -/*FUNCTION Icefront::CreatePVectorDiagnosticHO{{{*/ -ElementVector* Icefront::CreatePVectorDiagnosticHO(void){ - - /*Constants*/ - const int numdofs = NUMVERTICESQUA *NDOF2; - - /*Intermediaries*/ - int i,j,ig,index1,index2,index3,index4; - int fill; - IssmDouble surface,pressure,ice_pressure,rho_water,rho_ice,gravity; - IssmDouble water_pressure,air_pressure; - IssmDouble Jdet,z_g; - IssmDouble xyz_list[NUMVERTICESQUA][3]; - IssmDouble normal[3]; - IssmDouble l1l4[4]; - GaussPenta *gauss = NULL; - - Penta* penta=(Penta*)element; - - /*Initialize Element vector and return if necessary*/ - if(penta->IsOnWater()) return NULL; - ElementVector* pe=new ElementVector(nodes,NUMVERTICESQUA,this->parameters,HOApproximationEnum); - - /*Retrieve all inputs and parameters*/ - GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICESQUA); - Input* surface_input =penta->inputs->GetInput(SurfaceEnum); _assert_(surface_input); - inputs->GetInputValue(&fill,FillEnum); - rho_water=matpar->GetRhoWater(); - rho_ice =matpar->GetRhoIce(); - gravity =matpar->GetG(); - GetQuadNormal(&normal[0],xyz_list); - - /*Identify which nodes are in the quad: */ - index1=element->GetNodeIndex(nodes[0]); - index2=element->GetNodeIndex(nodes[1]); - index3=element->GetNodeIndex(nodes[2]); - index4=element->GetNodeIndex(nodes[3]); - - /* Start looping on the number of gaussian points: */ - IssmDouble zmax=xyz_list[0][2]; for(i=1;izmax) zmax=xyz_list[i][2]; - IssmDouble zmin=xyz_list[0][2]; for(i=1;i0 && zmin<0) gauss=new GaussPenta(index1,index2,index3,index4,3,10); //refined in vertical because of the sea level discontinuity - else gauss=new GaussPenta(index1,index2,index3,index4,3,3); - for(ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - penta->GetQuadNodalFunctions(l1l4,gauss,index1,index2,index3,index4); - penta->GetQuadJacobianDeterminant(&Jdet,xyz_list,gauss); - z_g=penta->GetZcoord(gauss); - surface_input->GetInputValue(&surface,gauss); - - switch(fill){ - case WaterEnum: - water_pressure=rho_water*gravity*min(0.,z_g);//0 if the gaussian point is above water level - break; - case AirEnum: - water_pressure=0; - break; - default: - _error_("fill type " << EnumToStringx(fill) << " not supported yet"); - } - ice_pressure=rho_ice*gravity*(surface-z_g); - air_pressure=0; - pressure = ice_pressure + water_pressure + air_pressure; - - for(i=0;ivalues[i*NDOF2+j]+=Jdet*gauss->weight*pressure*l1l4[i]*normal[j]; - } - - /*Transform load vector*/ - TransformLoadVectorCoord(pe,nodes,NUMVERTICESQUA,XYEnum); - - /*Clean up and return*/ - delete gauss; - return pe; -} -/*}}}*/ -/*FUNCTION Icefront::CreatePVectorDiagnosticFS{{{*/ -ElementVector* Icefront::CreatePVectorDiagnosticFS(void){ - - /*Constants*/ - const int numdofs = NUMVERTICESQUA *NDOF4; - - /*Intermediaries*/ - int i,j,ig,index1,index2,index3,index4; - int fill; - IssmDouble pressure,rho_water,gravity; - IssmDouble water_pressure,air_pressure; - IssmDouble Jdet,z_g; - IssmDouble xyz_list[NUMVERTICESQUA][3]; - IssmDouble normal[3]; - IssmDouble l1l4[4]; - GaussPenta *gauss = NULL; - - Penta* penta=(Penta*)element; - - /*Initialize Element vector and return if necessary*/ - if(penta->IsOnWater()) return NULL; - ElementVector* pe=new ElementVector(nodes,NUMVERTICESQUA,this->parameters,FSApproximationEnum); - - /*Retrieve all inputs and parameters*/ - GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICESQUA); - inputs->GetInputValue(&fill,FillEnum); - rho_water=matpar->GetRhoWater(); - gravity =matpar->GetG(); - GetQuadNormal(&normal[0],xyz_list); - - /*Identify which nodes are in the quad: */ - index1=element->GetNodeIndex(nodes[0]); - index2=element->GetNodeIndex(nodes[1]); - index3=element->GetNodeIndex(nodes[2]); - index4=element->GetNodeIndex(nodes[3]); - - /* Start looping on the number of gaussian points: */ - IssmDouble zmax=xyz_list[0][2]; for(i=1;izmax) zmax=xyz_list[i][2]; - IssmDouble zmin=xyz_list[0][2]; for(i=1;i0 && zmin<0) gauss=new GaussPenta(index1,index2,index3,index4,3,30); //refined in vertical because of the sea level discontinuity - else gauss=new GaussPenta(index1,index2,index3,index4,3,3); - for(ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - penta->GetQuadNodalFunctions(l1l4,gauss,index1,index2,index3,index4); - penta->GetQuadJacobianDeterminant(&Jdet,xyz_list,gauss); - z_g=penta->GetZcoord(gauss); - - switch(fill){ - case WaterEnum: - water_pressure=rho_water*gravity*min(0.,z_g);//0 if the gaussian point is above water level - break; - case AirEnum: - water_pressure=0; - break; - default: - _error_("fill type " << EnumToStringx(fill) << " not supported yet"); - } - air_pressure=0; - pressure = water_pressure + air_pressure; //no ice pressure fore FS - - for(i=0;ivalues[i*NDOF4+j]+=Jdet*gauss->weight*pressure*l1l4[i]*normal[j]; - else pe->values[i*NDOF4+j]+=0; //pressure term - } - } - } - - /*Transform load vector*/ - TransformLoadVectorCoord(pe,nodes,NUMVERTICESQUA,XYZPEnum); - - /*Clean up and return*/ - delete gauss; - return pe; -} -/*}}}*/ -#endif -/*FUNCTION Icefront::GetDofList {{{*/ -void Icefront::GetDofList(int** pdoflist,int approximation_enum,int setenum){ - - int numberofdofs=0; - int count=0; - int type; - int numberofnodes=2; - - /*output: */ - int* doflist=NULL; - - /*recover type: */ - inputs->GetInputValue(&type,IceFrontTypeEnum); - - /*Some checks for debugging*/ - _assert_(nodes); - - /*How many nodes? :*/ - if(type==SSA2dIceFrontEnum || type==SSA3dIceFrontEnum) - numberofnodes=2; - else - numberofnodes=4; - - /*Figure out size of doflist: */ - for(int i=0;iGetNumberOfDofs(approximation_enum,setenum); - } - - /*Allocate: */ - doflist=xNew(numberofdofs); - - /*Populate: */ - count=0; - for(int i=0;iGetDofList(doflist+count,approximation_enum,setenum); - count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum); - } - - /*Assign output pointers:*/ - *pdoflist=doflist; -} -/*}}}*/ -/*FUNCTION Icefront::GetSegmentNormal {{{*/ -void Icefront:: GetSegmentNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]){ - - /*Build unit outward pointing vector*/ - IssmDouble vector[2]; - IssmDouble norm; - - vector[0]=xyz_list[1][0] - xyz_list[0][0]; - vector[1]=xyz_list[1][1] - xyz_list[0][1]; - - norm=sqrt(pow(vector[0],2.0)+pow(vector[1],2.0)); - - normal[0]= + vector[1]/norm; - normal[1]= - vector[0]/norm; -} -/*}}}*/ -/*FUNCTION Icefront::GetQuadNormal {{{*/ -void Icefront:: GetQuadNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]){ - - /*Build unit outward pointing vector*/ - IssmDouble AB[3]; - IssmDouble AC[3]; - IssmDouble norm; - - AB[0]=xyz_list[1][0] - xyz_list[0][0]; - AB[1]=xyz_list[1][1] - xyz_list[0][1]; - AB[2]=xyz_list[1][2] - xyz_list[0][2]; - AC[0]=xyz_list[2][0] - xyz_list[0][0]; - AC[1]=xyz_list[2][1] - xyz_list[0][1]; - AC[2]=xyz_list[2][2] - xyz_list[0][2]; - - cross(normal,AB,AC); - norm=sqrt(pow(normal[0],2.0)+pow(normal[1],2.0)+pow(normal[2],2.0)); - - for(int i=0;i<3;i++) normal[i]=normal[i]/norm; -} -/*}}}*/ Index: sm/trunk-jpl/src/c/classes/Loads/Icefront.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Icefront.h (revision 15648) +++ (revision ) @@ -1,105 +1,0 @@ -/*!\file Icefront.h - * \brief: header file for icefront object - */ - -#ifndef _ICEFRONT_H_ -#define _ICEFRONT_H_ - -/*Headers:*/ -/*{{{*/ -#include "./Load.h" -class Hook; -class Inputs; -class Parameters; -class Matpar; -class Node; -class Element; -class IoModel; -class ElementVector; -class Vertex; -/*}}}*/ - -class Icefront: public Load { - - public: - int id; - int analysis_type; - - /*hooks: */ - Hook *hnodes; - Hook *hvertices; - Hook *helement; - Hook *hmatpar; - - /*Corresponding fields*/ - Matpar *matpar; - Node **nodes; - Vertex **vertices; - Element *element; - - /*inputs and parameters: */ - Inputs *inputs; - Parameters *parameters; - - /*Icefront constructors, destructors: {{{*/ - Icefront(); - Icefront(int icefront_id,int i, IoModel* iomodel,int in_icefront_type, int analysis_type); - Icefront(int element_id,const char* element_type_in,Inputs* inputs_in,Matpar* matpar_in, int icefront_type, int in_analysis_type); - ~Icefront(); - /*}}}*/ - /*Object virtual functions definitions:{{{ */ - void Echo(); - void DeepEcho(); - int Id(); - int ObjectEnum(); - Object* copy(); - /*}}}*/ - /*Update virtual functions definitions: {{{*/ - void InputUpdateFromVector(IssmDouble* vector, int name, int type); - void InputUpdateFromVector(int* vector, int name, int type); - void InputUpdateFromVector(bool* vector, int name, int type); - void InputUpdateFromMatrixDakota(IssmDouble* matrix,int ncols,int nrows, int name, int type); - void InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type); - void InputUpdateFromVectorDakota(int* vector, int name, int type); - void InputUpdateFromVectorDakota(bool* vector, int name, int type); - void InputUpdateFromConstant(IssmDouble constant, int name); - void InputUpdateFromConstant(int constant, int name); - void InputUpdateFromConstant(bool constant, int name); - void InputUpdateFromSolution(IssmDouble* solution); - void InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");}; - /*}}}*/ - /*Load virtual functions definitions: {{{*/ - void Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters); - void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters); - void CreateKMatrix(Matrix* Kff, Matrix* Kfs); - void CreatePVector(Vector* pf); - void CreateJacobianMatrix(Matrix* Jff); - int GetNumberOfNodes(void); - void GetNodesSidList(int* sidlist); - bool IsPenalty(void); - void PenaltyCreateKMatrix(Matrix* Kff, Matrix* kfs, IssmDouble kmax); - void PenaltyCreatePVector(Vector* pf, IssmDouble kmax); - void PenaltyCreateJacobianMatrix(Matrix* Jff,IssmDouble kmax); - void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int set1_enum,int set2_enum); - bool InAnalysis(int analysis_type); - /*}}}*/ - /*Load management: {{{*/ - void GetDofList(int** pdoflist,int approximation_enum,int setenum); - void GetSegmentNormal(IssmDouble* normal,IssmDouble xyz_list[2][3]); - void GetQuadNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]); - #ifdef _HAVE_CONTROL_ - ElementVector* CreatePVectorAdjointHoriz(void); - #endif - #ifdef _HAVE_DIAGNOSTIC_ - ElementVector* CreatePVectorDiagnosticHoriz(void); - ElementVector* CreatePVectorDiagnosticSSA2d(void); - #endif - #ifdef _HAVE_3D_ - ElementVector* CreatePVectorDiagnosticSSA3d(void); - ElementVector* CreatePVectorDiagnosticHO(void); - ElementVector* CreatePVectorDiagnosticFS(void); - #endif - /*}}}*/ -}; - -#endif /* _ICEFRONT_H_ */ Index: /issm/trunk-jpl/src/c/classes/Loads/Load.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Load.h (revision 15648) +++ /issm/trunk-jpl/src/c/classes/Loads/Load.h (revision 15649) @@ -9,5 +9,4 @@ /*Headers:*/ -/*{{{*/ class Node; template class Matrix; @@ -19,9 +18,6 @@ class Materials; class Parameters; - #include "../../datastructures/datastructures.h" #include "../Update.h" - -/*}}}*/ class Load: public Object,public Update{ Index: /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.h (revision 15648) +++ /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.h (revision 15649) @@ -7,5 +7,4 @@ /*Headers:*/ -/*{{{*/ #include "./Load.h" class Hook; @@ -13,7 +12,8 @@ class Inputs; class IoModel; +class Element; +class Vertex; class ElementMatrix; class ElementVector; -/*}}}*/ class Numericalflux: public Load {