Index: /issm/trunk/src/c/DataSet/DataSet.cpp =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.cpp (revision 1627) +++ /issm/trunk/src/c/DataSet/DataSet.cpp (revision 1628) @@ -300,4 +300,10 @@ rgb->Demarshall(&marshalled_dataset); dataset->AddObject(rgb); + } + else if(enum_type==RiftfrontEnum()){ + Riftfront* riftfront=NULL; + riftfront=new Riftfront(); + riftfront->Demarshall(&marshalled_dataset); + dataset->AddObject(riftfront); } else{ @@ -1371,6 +1377,10 @@ } } - - - -} +} + + +int DataSet::GetEnum(int offset){ + + return objects[offset]->Enum(); + +} Index: /issm/trunk/src/c/DataSet/DataSet.h =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.h (revision 1627) +++ /issm/trunk/src/c/DataSet/DataSet.h (revision 1628) @@ -36,4 +36,5 @@ int GetEnum(); + int GetEnum(int offset); void Echo(); void DeepEcho(); Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp (revision 1627) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp (revision 1628) @@ -70,8 +70,10 @@ int NodeEnum(void){ return 420; } /*Loads: */ -int LoadEnum(void){ return 430; } -int IcefrontEnum(void){ return 431; } -int PenpairEnum(void){ return 432; } -int PengridEnum(void){ return 433; } +int LoadEnum(void){ return 200; } +int IcefrontEnum(void){ return 201; } +int RiftfrontEnum(void){ return 202; } +int PenpairEnum(void){ return 203; } +int PengridEnum(void){ return 204; } + /*Materials: */ int MaterialEnum(void){ return 440; } Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 1627) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 1628) @@ -20,4 +20,5 @@ int PengridEnum(void); int IcefrontEnum(void); +int RiftfrontEnum(void); int ParamEnum(void); int ElementEnum(void); Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 1627) +++ /issm/trunk/src/c/Makefile.am (revision 1628) @@ -64,4 +64,6 @@ ./objects/Icefront.cpp\ ./objects/Icefront.h\ + ./objects/Riftfront.cpp\ + ./objects/Riftfront.h\ ./objects/Param.cpp\ ./objects/Param.h\ @@ -169,5 +171,4 @@ ./io/WriteParams.cpp\ ./io/FetchNodeSets.cpp\ - ./io/FetchRifts.cpp\ ./io/ParameterInputsInit.cpp\ ./io/pfopen.cpp\ @@ -272,4 +273,6 @@ ./PenaltyConstraintsx/PenaltyConstraintsx.cpp\ ./PenaltyConstraintsx/PenaltyConstraintsx.h\ + ./PenaltyConstraintsx/RiftConstraints.cpp\ + ./PenaltyConstraintsx/RiftConstraints.h\ ./PenaltySystemMatricesx/PenaltySystemMatricesx.cpp\ ./PenaltySystemMatricesx/PenaltySystemMatricesx.h\ @@ -350,4 +353,6 @@ ./objects/Icefront.cpp\ ./objects/Icefront.h\ + ./objects/Riftfront.cpp\ + ./objects/Riftfront.h\ ./objects/Param.cpp\ ./objects/Param.h\ @@ -453,5 +458,4 @@ ./io/WriteParams.cpp\ ./io/FetchNodeSets.cpp\ - ./io/FetchRifts.cpp\ ./io/ParameterInputsInit.cpp\ ./io/pfopen.cpp\ @@ -550,4 +554,6 @@ ./PenaltyConstraintsx/PenaltyConstraintsx.cpp\ ./PenaltyConstraintsx/PenaltyConstraintsx.h\ + ./PenaltyConstraintsx/RiftConstraints.cpp\ + ./PenaltyConstraintsx/RiftConstraints.h\ ./PenaltySystemMatricesx/PenaltySystemMatricesx.cpp\ ./PenaltySystemMatricesx/PenaltySystemMatricesx.h\ Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp (revision 1627) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp (revision 1628) @@ -143,9 +143,4 @@ /*Rifts:*/ - int* riftsnumpenaltypairs; - double** riftspenaltypairs; - int* riftsfill; - double* riftsfriction; - double* riftpenaltypairs=NULL; int el1,el2; @@ -193,23 +188,13 @@ /*Deal with rifts, they have to be included into one partition only, not several: */ - FetchRifts(&riftsnumpenaltypairs,&riftspenaltypairs,&riftsfill,&riftsfriction,model_handle,model->numrifts); + ModelFetchData((void**)&model->riftinfo,&model->numrifts,NULL,model_handle,"riftinfo","Matrix","Mat"); for(i=0;inumrifts;i++){ - riftpenaltypairs=model->riftspenaltypairs[i]; - for(j=0;jriftsnumpenaltypairs[i];j++){ - el1=(int)*(riftpenaltypairs+7*j+2)-1; //matlab indexing to c indexing - el2=(int)*(riftpenaltypairs+7*j+3)-1; //matlab indexing to c indexing - epart[el2]=epart[el1]; //ensures that this pair of elements will be in the same partition, as well as the corresponding grids; - } + el1=(int)*(model->riftinfo+9*j+2)-1; //matlab indexing to c indexing + el2=(int)*(model->riftinfo+9*j+3)-1; //matlab indexing to c indexing + epart[el2]=epart[el1]; //ensures that this pair of elements will be in the same partition, as well as the corresponding grids; } /*Free rifts: */ - xfree((void**)&riftsnumpenaltypairs); - for(i=0;inumrifts;i++){ - double* temp=riftspenaltypairs[i]; - xfree((void**)&temp); - } - xfree((void**)&riftspenaltypairs); - xfree((void**)&riftsfill); - xfree((void**)&riftsfriction); + xfree((void**)&model->riftinfo); /*Used later on: */ Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp (revision 1627) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp (revision 1628) @@ -16,5 +16,5 @@ void CreateLoadsDiagnosticHoriz(DataSet** ploads, Model* model,ConstDataHandle model_handle){ - int i,j,counter; + int i; int element; @@ -24,5 +24,5 @@ DataSet* loads = NULL; Icefront* icefront = NULL; - Penpair* penpair = NULL; + Riftfront* riftfront= NULL; int segment_width; @@ -30,11 +30,5 @@ int i1,i2,i3,i4; int i0,range; - int grid1,grid2; - - /*rift penpair: */ - double normal[2]; - double length; - double friction; - int fill; + /*icefront intermediary data: */ @@ -48,24 +42,28 @@ double icefront_b[MAX_ICEFRONT_GRIDS]; - /*penpair intermediary data: */ - int penpair_id; - int penpair_node_ids[2]; - double penpair_penalty_offset; - int penpair_numdofs; - int penpair_dof; - int penpair_penalty_lock; - int penpair_element_ids[2]; - double penpair_friction; - int penpair_fill; - double penpair_normal[2]; - double penpair_length; - - /*Rifts:*/ - int* riftsnumpenaltypairs=NULL; - double** riftspenaltypairs=NULL; - int* riftsfill=NULL; - int* riftsfriction=NULL; - double* riftpenaltypairs=NULL; + /*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; + bool riftfront_shelf; + double riftfront_penalty_offset; + int riftfront_penalty_lock; + bool riftfront_active; + int riftfront_counter; + bool riftfront_prestable; int el1,el2; + int grid1,grid2; + double normal[2]; + double length; + int fill; + double friction; /*Create loads: */ @@ -153,56 +151,81 @@ xfree((void**)&model->bed); - counter=0; - /*Create penpair loads also for rift grids, so that they don't penetrate one another, if needed: */ /*First fetch rifts: */ + ModelFetchData((void**)&model->riftinfo,&model->numrifts,NULL,model_handle,"riftinfo","Matrix","Mat"); + ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat"); + ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat"); + ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat"); + ModelFetchData((void**)&model->gridoniceshelf,NULL,NULL,model_handle,"gridoniceshelf","Matrix","Mat"); + if(model->numrifts){ - for(i=0;inumrifts;i++){ - riftpenaltypairs=model->riftspenaltypairs[i]; - for(j=0;jriftsnumpenaltypairs[i];j++){ - - el1=(int)*(riftpenaltypairs+7*j+2); - #ifdef _PARALLEL_ - if (model->epart[el1-1]!=my_rank){ - /*This load does not belong to this cluster node, as it references an element - *that does not belong to this node's partition. Drop this 'j':*/ - continue; - } - #endif - - /*Ok, retrieve all the data needed to add a penalty between the two grids: */ - el2=(int)*(riftpenaltypairs+7*j+3); - grid1=(int)*(riftpenaltypairs+7*j+0); - grid2=(int)*(riftpenaltypairs+7*j+1); - normal[0]=*(riftpenaltypairs+7*j+4); - normal[1]=*(riftpenaltypairs+7*j+5); - length=*(riftpenaltypairs+7*j+6); - friction=model->riftsfriction[i]; - fill=model->riftsfill[i]; - - penpair_id=counter+1; //matlab indexing - penpair_node_ids[0]=grid1; - penpair_node_ids[1]=grid2; - penpair_element_ids[0]=el1; - penpair_element_ids[1]=el2; - penpair_penalty_offset=model->penalty_offset; - penpair_penalty_lock=model->penalty_lock; - penpair_numdofs=2; - penpair_normal[0]=normal[0]; - penpair_normal[1]=normal[1]; - penpair_length=length; - penpair_friction=friction; - penpair_fill=fill; - - penpair= new Penpair(penpair_id,penpair_penalty_offset,penpair_penalty_lock,penpair_numdofs,penpair_node_ids,penpair_dof, - penpair_element_ids,penpair_friction,penpair_fill,penpair_normal,penpair_length); - loads->AddObject(penpair); - - counter++; + el1=(int)*(model->riftinfo+9*i+2); + #ifdef _PARALLEL_ + if (model->epart[el1-1]!=my_rank){ + /*This load does not belong to this cluster node, as it references an element + *that does not belong to this node's partition. Drop this 'j':*/ + continue; } + #endif + + /*Ok, retrieve all the data needed to add a penalty between the two grids: */ + el2=(int)*(model->riftinfo+9*i+3); + + grid1=(int)*(model->riftinfo+9*i+0); + grid2=(int)*(model->riftinfo+9*i+1); + + normal[0]=*(model->riftinfo+9*i+4); + normal[1]=*(model->riftinfo+9*i+5); + length=*(model->riftinfo+9*i+6); + + fill = (int)*(model->riftinfo+9*i+7); + friction=*(model->riftinfo+9*i+8); + + strcpy(riftfront_type,"2d"); + riftfront_id=i+1; //matlab indexing + riftfront_node_ids[0]=grid1; + riftfront_node_ids[1]=grid2; + riftfront_mparid=model->numberofelements+1; //matlab indexing + + riftfront_h[0]=model->thickness[grid1-1]; + riftfront_h[1]=model->thickness[grid2-1]; + + riftfront_b[0]=model->bed[grid1-1]; + riftfront_b[1]=model->bed[grid2-1]; + + riftfront_s[0]=model->surface[grid1-1]; + riftfront_s[1]=model->surface[grid2-1]; + + riftfront_normal[0]=normal[0]; + riftfront_normal[1]=normal[1]; + riftfront_length=length; + + riftfront_fill=fill; + riftfront_friction=friction; + riftfront_shelf=(bool)model->gridoniceshelf[grid1-1]; + + riftfront_penalty_offset=model->penalty_offset; + riftfront_penalty_lock=model->penalty_lock; + + riftfront_active=0; + riftfront_counter=0; + riftfront_prestable=0; + + riftfront=new Riftfront(riftfront_type,riftfront_id, riftfront_node_ids, riftfront_mparid, riftfront_h,riftfront_b,riftfront_s,riftfront_normal,riftfront_length,riftfront_fill,riftfront_friction, riftfront_penalty_offset, riftfront_penalty_lock, riftfront_active,riftfront_counter,riftfront_prestable,riftfront_shelf); + + loads->AddObject(riftfront); } } + + /*free ressources: */ + xfree((void**)&model->riftinfo); + xfree((void**)&model->thickness); + xfree((void**)&model->bed); + xfree((void**)&model->surface); + xfree((void**)&model->gridoniceshelf); + + cleanup_and_return: /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these @@ -210,17 +233,4 @@ loads->Presort(); - - /*Free ressources:*/ - xfree((void**)&riftsnumpenaltypairs); - for(i=0;inumrifts;i++){ - double* temp=riftspenaltypairs[i]; - xfree((void**)&temp); - } - xfree((void**)&riftspenaltypairs); - xfree((void**)&riftsfill); - xfree((void**)&riftsfriction); - - cleanup_and_return: - /*Assign output pointer: */ *ploads=loads; Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp (revision 1627) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp (revision 1628) @@ -121,12 +121,4 @@ int first_grid_index; - /*Rifts:*/ - int* riftsnumpenaltypairs; - double** riftspenaltypairs; - int* riftsfill; - double* riftsfriction; - double* riftpenaltypairs=NULL; - int el1,el2; - /*First create the elements, nodes and material properties: */ elements = new DataSet(ElementsEnum()); @@ -167,26 +159,4 @@ xfree((void**)&model->elements); xfree((void**)&model->elements2d); - - - /*Deal with rifts, they have to be included into one partition only, not several: */ - FetchRifts(&riftsnumpenaltypairs,&riftspenaltypairs,&riftsfill,&riftsfriction,model_handle,model->numrifts); - - for(i=0;inumrifts;i++){ - riftpenaltypairs=model->riftspenaltypairs[i]; - for(j=0;jriftsnumpenaltypairs[i];j++){ - el1=(int)*(riftpenaltypairs+7*j+2)-1; //matlab indexing to c indexing - el2=(int)*(riftpenaltypairs+7*j+3)-1; //matlab indexing to c indexing - epart[el2]=epart[el1]; //ensures that this pair of elements will be in the same partition, as well as the corresponding grids; - } - } - /*Free rifts: */ - xfree((void**)&riftsnumpenaltypairs); - for(i=0;inumrifts;i++){ - double* temp=riftspenaltypairs[i]; - xfree((void**)&temp); - } - xfree((void**)&riftspenaltypairs); - xfree((void**)&riftsfill); - xfree((void**)&riftsfriction); /*Used later on: */ Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp (revision 1627) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp (revision 1628) @@ -120,12 +120,4 @@ int first_grid_index; - /*Rifts:*/ - int* riftsnumpenaltypairs; - double** riftspenaltypairs; - int* riftsfill; - double* riftsfriction; - double* riftpenaltypairs=NULL; - int el1,el2; - /*Penalty partitioning: */ int num_grids3d_collapsed; @@ -167,26 +159,4 @@ xfree((void**)&model->elements2d); - - /*Deal with rifts, they have to be included into one partition only, not several: */ - FetchRifts(&riftsnumpenaltypairs,&riftspenaltypairs,&riftsfill,&riftsfriction,model_handle,model->numrifts); - - for(i=0;inumrifts;i++){ - riftpenaltypairs=model->riftspenaltypairs[i]; - for(j=0;jriftsnumpenaltypairs[i];j++){ - el1=(int)*(riftpenaltypairs+7*j+2)-1; //matlab indexing to c indexing - el2=(int)*(riftpenaltypairs+7*j+3)-1; //matlab indexing to c indexing - epart[el2]=epart[el1]; //ensures that this pair of elements will be in the same partition, as well as the corresponding grids; - } - } - /*Free rifts: */ - xfree((void**)&riftsnumpenaltypairs); - for(i=0;inumrifts;i++){ - double* temp=riftspenaltypairs[i]; - xfree((void**)&temp); - } - xfree((void**)&riftspenaltypairs); - xfree((void**)&riftsfill); - xfree((void**)&riftsfriction); - /*Used later on: */ my_grids=(int*)xcalloc(model->numberofnodes,sizeof(int)); Index: /issm/trunk/src/c/ModelProcessorx/Melting/CreateElementsNodesAndMaterialsMelting.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Melting/CreateElementsNodesAndMaterialsMelting.cpp (revision 1627) +++ /issm/trunk/src/c/ModelProcessorx/Melting/CreateElementsNodesAndMaterialsMelting.cpp (revision 1628) @@ -120,12 +120,4 @@ int first_grid_index; - /*Rifts:*/ - int* riftsnumpenaltypairs; - double** riftspenaltypairs; - int* riftsfill; - double* riftsfriction; - double* riftpenaltypairs=NULL; - int el1,el2; - /*First create the elements, nodes and material properties: */ elements = new DataSet(ElementsEnum()); @@ -151,26 +143,4 @@ /*Free elements and elements2d: */ xfree((void**)&model->elements2d); - - - /*Deal with rifts, they have to be included into one partition only, not several: */ - FetchRifts(&riftsnumpenaltypairs,&riftspenaltypairs,&riftsfill,&riftsfriction,model_handle,model->numrifts); - - for(i=0;inumrifts;i++){ - riftpenaltypairs=model->riftspenaltypairs[i]; - for(j=0;jriftsnumpenaltypairs[i];j++){ - el1=(int)*(riftpenaltypairs+7*j+2)-1; //matlab indexing to c indexing - el2=(int)*(riftpenaltypairs+7*j+3)-1; //matlab indexing to c indexing - epart[el2]=epart[el1]; //ensures that this pair of elements will be in the same partition, as well as the corresponding grids; - } - } - /*Free rifts: */ - xfree((void**)&riftsnumpenaltypairs); - for(i=0;inumrifts;i++){ - double* temp=riftspenaltypairs[i]; - xfree((void**)&temp); - } - xfree((void**)&riftspenaltypairs); - xfree((void**)&riftsfill); - xfree((void**)&riftsfriction); /*Used later on: */ Index: /issm/trunk/src/c/ModelProcessorx/Model.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Model.cpp (revision 1627) +++ /issm/trunk/src/c/ModelProcessorx/Model.cpp (revision 1628) @@ -156,8 +156,5 @@ model->numrifts=0; - model->riftsnumpenaltypairs=NULL; - model->riftspenaltypairs=NULL; - model->riftsfill=NULL; - model->riftsfriction=NULL; + model->riftinfo=NULL; /*!penalties: */ @@ -271,14 +268,5 @@ xfree((void**)&model->sub_analysis_type); - if(model->numrifts){ - for(i=0;inumrifts;i++){ - double* riftpenaltypairs=model->riftspenaltypairs[i]; - xfree((void**)&riftpenaltypairs); - } - xfree((void**)&model->riftspenaltypairs); - xfree((void**)&model->riftsnumpenaltypairs); - xfree((void**)&model->riftsfill); - xfree((void**)&model->riftsfriction); - } + xfree((void**)&model->riftinfo); xfree((void**)&model->penalties); @@ -387,7 +375,4 @@ ModelFetchData((void**)&model->thermal_exchange_velocity,NULL,NULL,model_handle,"thermal_exchange_velocity","Scalar",NULL); - /*rifts: */ - ModelFetchData((void**)&model->numrifts,NULL,NULL,model_handle,"numrifts","Integer",NULL); - /*qmu: */ if(model->qmu_analysis){ @@ -426,18 +411,4 @@ } - if(which_part==2 && my_rank==rank){ - printf("Model rifts: \n"); - printf(" number of rifts: %i\n",model->numrifts); - for(i=0;inumrifts;i++){ - double* penaltypairs=model->riftspenaltypairs[i]; - printf(" rift #%i\n",i); - for (j=0;jriftsnumpenaltypairs[i];j++){ - printf(" grids %g %g elements %g %g normal [%g,%g] length %g\n",*(penaltypairs+7*j+0),*(penaltypairs+7*j+1),*(penaltypairs+7*j+2),*(penaltypairs+7*j+3), - *(penaltypairs+7*j+4),*(penaltypairs+7*j+5),*(penaltypairs+7*j+6)); - } - printf(" friction %g fill %i\n",model->riftsfriction[i],model->riftsfill[i]); - } - } - cleanup_and_return: return; } Index: /issm/trunk/src/c/ModelProcessorx/Model.h =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Model.h (revision 1627) +++ /issm/trunk/src/c/ModelProcessorx/Model.h (revision 1628) @@ -152,8 +152,5 @@ /*rifts: */ int numrifts; - int* riftsnumpenaltypairs; - double** riftspenaltypairs; - int* riftsfill; - double* riftsfriction; + double* riftinfo; /*penalties: */ Index: /issm/trunk/src/c/ModelProcessorx/Prognostic/CreateElementsNodesAndMaterialsPrognostic.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Prognostic/CreateElementsNodesAndMaterialsPrognostic.cpp (revision 1627) +++ /issm/trunk/src/c/ModelProcessorx/Prognostic/CreateElementsNodesAndMaterialsPrognostic.cpp (revision 1628) @@ -143,12 +143,4 @@ int first_grid_index; - /*Rifts:*/ - int* riftsnumpenaltypairs; - double** riftspenaltypairs; - int* riftsfill; - double* riftsfriction; - double* riftpenaltypairs=NULL; - int el1,el2; - /*First create the elements, nodes and material properties: */ elements = new DataSet(ElementsEnum()); @@ -187,26 +179,4 @@ xfree((void**)&model->elements); xfree((void**)&model->elements2d); - - - /*Deal with rifts, they have to be included into one partition only, not several: */ - FetchRifts(&riftsnumpenaltypairs,&riftspenaltypairs,&riftsfill,&riftsfriction,model_handle,model->numrifts); - - for(i=0;inumrifts;i++){ - riftpenaltypairs=model->riftspenaltypairs[i]; - for(j=0;jriftsnumpenaltypairs[i];j++){ - el1=(int)*(riftpenaltypairs+7*j+2)-1; //matlab indexing to c indexing - el2=(int)*(riftpenaltypairs+7*j+3)-1; //matlab indexing to c indexing - epart[el2]=epart[el1]; //ensures that this pair of elements will be in the same partition, as well as the corresponding grids; - } - } - /*Free rifts: */ - xfree((void**)&riftsnumpenaltypairs); - for(i=0;inumrifts;i++){ - double* temp=riftspenaltypairs[i]; - xfree((void**)&temp); - } - xfree((void**)&riftspenaltypairs); - xfree((void**)&riftsfill); - xfree((void**)&riftsfriction); /*Used later on: */ Index: /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp (revision 1627) +++ /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp (revision 1628) @@ -121,12 +121,4 @@ int first_grid_index; - /*Rifts:*/ - int* riftsnumpenaltypairs; - double** riftspenaltypairs; - int* riftsfill; - double* riftsfriction; - double* riftpenaltypairs=NULL; - int el1,el2; - /*First create the elements, nodes and material properties: */ elements = new DataSet(ElementsEnum()); @@ -167,25 +159,4 @@ xfree((void**)&model->elements2d); - - /*Deal with rifts, they have to be included into one partition only, not several: */ - FetchRifts(&riftsnumpenaltypairs,&riftspenaltypairs,&riftsfill,&riftsfriction,model_handle,model->numrifts); - - for(i=0;inumrifts;i++){ - riftpenaltypairs=model->riftspenaltypairs[i]; - for(j=0;jriftsnumpenaltypairs[i];j++){ - el1=(int)*(riftpenaltypairs+7*j+2)-1; //matlab indexing to c indexing - el2=(int)*(riftpenaltypairs+7*j+3)-1; //matlab indexing to c indexing - epart[el2]=epart[el1]; //ensures that this pair of elements will be in the same partition, as well as the corresponding grids; - } - } - /*Free rifts: */ - xfree((void**)&riftsnumpenaltypairs); - for(i=0;inumrifts;i++){ - double* temp=riftspenaltypairs[i]; - xfree((void**)&temp); - } - xfree((void**)&riftspenaltypairs); - xfree((void**)&riftsfill); - xfree((void**)&riftsfriction); /*Used later on: */ Index: /issm/trunk/src/c/ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp (revision 1627) +++ /issm/trunk/src/c/ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp (revision 1628) @@ -121,12 +121,4 @@ int first_grid_index; - /*Rifts:*/ - int* riftsnumpenaltypairs; - double** riftspenaltypairs; - int* riftsfill; - double* riftsfriction; - double* riftpenaltypairs=NULL; - int el1,el2; - /*First create the elements, nodes and material properties: */ elements = new DataSet(ElementsEnum()); @@ -154,25 +146,4 @@ /*Free elements and elements2d: */ xfree((void**)&model->elements2d); - - /*Deal with rifts, they have to be included into one partition only, not several: */ - FetchRifts(&riftsnumpenaltypairs,&riftspenaltypairs,&riftsfill,&riftsfriction,model_handle,model->numrifts); - - for(i=0;inumrifts;i++){ - riftpenaltypairs=model->riftspenaltypairs[i]; - for(j=0;jriftsnumpenaltypairs[i];j++){ - el1=(int)*(riftpenaltypairs+7*j+2)-1; //matlab indexing to c indexing - el2=(int)*(riftpenaltypairs+7*j+3)-1; //matlab indexing to c indexing - epart[el2]=epart[el1]; //ensures that this pair of elements will be in the same partition, as well as the corresponding grids; - } - } - /*Free rifts: */ - xfree((void**)&riftsnumpenaltypairs); - for(i=0;inumrifts;i++){ - double* temp=riftspenaltypairs[i]; - xfree((void**)&temp); - } - xfree((void**)&riftspenaltypairs); - xfree((void**)&riftsfill); - xfree((void**)&riftsfriction); /*Used later on: */ Index: /issm/trunk/src/c/PenaltyConstraintsx/PenaltyConstraintsx.cpp =================================================================== --- /issm/trunk/src/c/PenaltyConstraintsx/PenaltyConstraintsx.cpp (revision 1627) +++ /issm/trunk/src/c/PenaltyConstraintsx/PenaltyConstraintsx.cpp (revision 1628) @@ -4,4 +4,5 @@ #include "./PenaltyConstraintsx.h" +#include "./RiftConstraints.h" #undef __FUNCT__ @@ -30,6 +31,6 @@ /*Do we have penalties linked to rifts? In this case, run our special rifts penalty * management routine, otherwise, skip : */ - if (loads->RiftIsPresent()){ - loads->RiftConstraints(&num_unstable_constraints,inputs,analysis_type,sub_analysis_type); + if (RiftIsPresent(loads)){ + RiftConstraints(&num_unstable_constraints,loads,inputs,analysis_type,sub_analysis_type); } else if(loads->MeltingIsPresent()){ Index: /issm/trunk/src/c/PenaltyConstraintsx/RiftConstraints.cpp =================================================================== --- /issm/trunk/src/c/PenaltyConstraintsx/RiftConstraints.cpp (revision 1628) +++ /issm/trunk/src/c/PenaltyConstraintsx/RiftConstraints.cpp (revision 1628) @@ -0,0 +1,297 @@ +/*!\file RiftConstraints.cpp + * \brief: manage penalties for rifts + */ + +#include "./RiftConstraints.h" +#include "../EnumDefinitions/EnumDefinitions.h" + +#undef __FUNCT__ +#define __FUNCT__ "RiftConstrain" +int RiftConstraints(int* pnum_unstable_constraints,DataSet* loads,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){ + + int num_unstable_constraints=0; + int potential; + + /*First, is set of constraints pre-stabilised? Ie: all unstable rifts are penalised? :*/ + if(!IsPreStable(loads)){ + + PreConstrain(&num_unstable_constraints,loads,inputs,analysis_type); + + if(!num_unstable_constraints){ + /*Ok, preconstraining is done. Figure out size of unstable constraints pool. Then set + * prestable flag to 1 for all rift grids: */ + + potential=PotentialUnstableConstraints(loads,inputs,analysis_type); + printf(" set of constraints is prestabilised; unstable constraints potentialnumber: %i\n",potential); + + SetPreStable(loads); + num_unstable_constraints=-1; //ensure convergence does not happen! + } + + } + else{ + /*Ok, set of constraints is stable. Start increasing number of inactive constraints, one by one: */ + + /*Figure out max penetration: */ + MaxPenetrationInInputs(loads,inputs,analysis_type); + + /*Constrain rifts: */ + Constrain(&num_unstable_constraints,loads,inputs,analysis_type); + } + + /*Assign output pointers: */ + *pnum_unstable_constraints=num_unstable_constraints; + +} + +#undef __FUNCT__ +#define __FUNCT__ "RiftIsPresent" + +int RiftIsPresent(DataSet* loads){ + + + int i; + + int found=0; + int mpi_found=0; + + /*go though loads, and figure out if one of the loads is a Riftfront: */ + for (i=0;iSize();i++){ + + if(RiftfrontEnum()==loads->GetEnum(i)){ + found=1; + break; + } + } + + #ifdef _PARALLEL_ + MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD ); + MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD); + found=mpi_found; + #endif + + return found; +} + +#undef __FUNCT__ +#define __FUNCT__ "IsPreStable" + +int IsPreStable(DataSet* loads){ + + + int i; + + Riftfront* riftfront=NULL; + int found=0; + int mpi_found=0; + + /*go though loads, and figure out if one of the penpair loads is still not stable: */ + for (i=0;iSize();i++){ + + if(RiftfrontEnum()==loads->GetEnum(i)){ + + riftfront=(Riftfront*)loads->GetObjectByOffset(i); + + if (riftfront->PreStable()==0){ + found=1; + break; + } + } + } + + #ifdef _PARALLEL_ + MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD ); + MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD); + found=mpi_found; + #endif + + if (found){ + /*We found an unstable constraint. : */ + return 0; + } + else{ + return 1; + } +} + +#undef __FUNCT__ +#define __FUNCT__ "SetPreStable" + +int SetPreStable(DataSet* loads){ + + + int i; + + Riftfront* riftfront=NULL; + int found=0; + int mpi_found=0; + + /*go though loads, and set loads to pre stable.:*/ + for (i=0;iSize();i++){ + + if(RiftfrontEnum()==loads->GetEnum(i)){ + + riftfront=(Riftfront*)loads->GetObjectByOffset(i); + riftfront->SetPreStable(); + } + } +} + +#undef __FUNCT__ +#define __FUNCT__ "PreConstrain" + +int PreConstrain(int* pnum_unstable_constraints,DataSet* loads,ParameterInputs* inputs,int analysis_type){ + + int i; + + /* generic object pointer: */ + Riftfront* riftfront=NULL; + + int unstable; + int sum_num_unstable_constraints; + int num_unstable_constraints=0; + + /*Enforce constraints: */ + for (i=0;iSize();i++){ + + if (RiftfrontEnum()==loads->GetEnum(i)){ + + riftfront=(Riftfront*)loads->GetObjectByOffset(i); + + riftfront->PreConstrain(&unstable,inputs,analysis_type); + + num_unstable_constraints+=unstable; + } + } + + #ifdef _PARALLEL_ + MPI_Reduce (&num_unstable_constraints,&sum_num_unstable_constraints,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD ); + MPI_Bcast(&sum_num_unstable_constraints,1,MPI_INT,0,MPI_COMM_WORLD); + num_unstable_constraints=sum_num_unstable_constraints; + #endif + + /*Assign output pointers: */ + *pnum_unstable_constraints=num_unstable_constraints; + +} + +#undef __FUNCT__ +#define __FUNCT__ "Constrain" + +int Constrain(int* pnum_unstable_constraints,DataSet* loads,ParameterInputs* inputs,int analysis_type){ + + int i; + + /* generic object pointer: */ + Riftfront* riftfront=NULL; + + int unstable; + int sum_num_unstable_constraints; + int num_unstable_constraints=0; + + /*Enforce constraints: */ + for (i=0;iSize();i++){ + + if (RiftfrontEnum()==loads->GetEnum(i)){ + + riftfront=(Riftfront*)loads->GetObjectByOffset(i); + + riftfront->Constrain(&unstable,inputs,analysis_type); + + num_unstable_constraints+=unstable; + } + } + + #ifdef _PARALLEL_ + MPI_Reduce (&num_unstable_constraints,&sum_num_unstable_constraints,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD ); + MPI_Bcast(&sum_num_unstable_constraints,1,MPI_INT,0,MPI_COMM_WORLD); + num_unstable_constraints=sum_num_unstable_constraints; + #endif + + /*Assign output pointers: */ + *pnum_unstable_constraints=num_unstable_constraints; + +} + +#undef __FUNCT__ +#define __FUNCT__ "MaxPenetrationInInputs" + +int MaxPenetrationInInputs(DataSet* loads,ParameterInputs* inputs,int analysis_type){ + + int i; + + /* generic object pointer: */ + Riftfront* riftfront=NULL; + + /*rift penetration: */ + double max_penetration=0; + double mpi_max_penetration; + double penetration; + + /*Ok, we are going to find the grid pairs which are not penetrating, even though they + * are penalised. We will release only the one with has least <0 penetration. : */ + + max_penetration=0; + for (i=0;iSize();i++){ + + if (RiftfrontEnum()==loads->GetEnum(i)){ + + riftfront=(Riftfront*)loads->GetObjectByOffset(i); + + riftfront->Penetration(&penetration,inputs,analysis_type); + + if (penetration>max_penetration)max_penetration=penetration; + } + } + + #ifdef _PARALLEL_ + MPI_Reduce (&max_penetration,&mpi_max_penetration,1,MPI_DOUBLE,MPI_MAX,0,MPI_COMM_WORLD ); + MPI_Bcast(&mpi_max_penetration,1,MPI_DOUBLE,0,MPI_COMM_WORLD); + max_penetration=mpi_max_penetration; + #endif + + #ifdef _DEBUG_ + printf("Maximum rift penetration2: %g\n",max_penetration); + #endif + + /*feed max_penetration to inputs: */ + inputs->Add("max_penetration",max_penetration); +} + + +#undef __FUNCT__ +#define __FUNCT__ "PotentialUnstableConstraints" + +int PotentialUnstableConstraints(DataSet* loads,ParameterInputs* inputs,int analysis_type){ + + int i; + + /* generic object pointer: */ + Riftfront* riftfront=NULL; + + /*Ok, we are going to find the grid pairs which are not penetrating, even though they + * are penalised. We will release only the one with has least <0 penetration. : */ + int unstable=0; + int sum_num_unstable_constraints=0; + int num_unstable_constraints=0; + + for (i=0;iSize();i++){ + + if (RiftfrontEnum()==loads->GetEnum(i)){ + + riftfront=(Riftfront*)loads->GetObjectByOffset(i); + + riftfront->PotentialUnstableConstraint(&unstable,inputs,analysis_type); + + num_unstable_constraints+=unstable; + } + } + + #ifdef _PARALLEL_ + MPI_Reduce (&num_unstable_constraints,&sum_num_unstable_constraints,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD ); + MPI_Bcast(&sum_num_unstable_constraints,1,MPI_INT,0,MPI_COMM_WORLD); + num_unstable_constraints=sum_num_unstable_constraints; + #endif + + return num_unstable_constraints; +} Index: /issm/trunk/src/c/PenaltyConstraintsx/RiftConstraints.h =================================================================== --- /issm/trunk/src/c/PenaltyConstraintsx/RiftConstraints.h (revision 1628) +++ /issm/trunk/src/c/PenaltyConstraintsx/RiftConstraints.h (revision 1628) @@ -0,0 +1,27 @@ +/*!\file RiftConstraints.h + * \brief: manage penalties for rifts + */ + + +#ifndef _RIFTCONSTRAINTS_H_ +#define _RIFTCONSTRAINTS_H_ + +#include "../objects/objects.h" +#include "../DataSet/DataSet.h" + +int RiftConstraints(int* pnum_unstable_constraints,DataSet* loads,ParameterInputs* inputs,int analysis_type,int sub_analysis_type); + +int RiftIsPresent(DataSet* loads); + +int IsPreStable(DataSet* loads); + +int SetPreStable(DataSet* loads); + +int PreConstrain(int* pnum_unstable_constraints,DataSet* loads,ParameterInputs* inputs,int analysis_type_enum); + +int Constrain(int* pnum_unstable_constraints,DataSet* loads,ParameterInputs* inputs,int analysis_type_enum); + +int MaxPenetrationInInputs(DataSet* loads,ParameterInputs* inputs,int analysis_type_enum); + +int PotentialUnstableConstraints(DataSet* loads,ParameterInputs* inputs,int analysis_type_enum); +#endif Index: sm/trunk/src/c/io/FetchRifts.cpp =================================================================== --- /issm/trunk/src/c/io/FetchRifts.cpp (revision 1627) +++ (revision ) @@ -1,78 +1,0 @@ -/*! \file FetchRifts.cpp - * \brief: special i/o here to recover the rifts data. Because serially, the rifts are held into a matlab array strucutre, - * hard to parse. In parallel, it's easier, as we have marshalled the data into a binary buffer. - */ - -#undef __FUNCT__ -#define __FUNCT__ "FetchRifts" - -#include "./io.h" -#include "../shared/shared.h" - -int FetchRifts(int** priftsnumpenaltypairs,double*** priftspenaltypairs,int** priftsfill,double** priftsfriction,ConstDataHandle model_handle,int numrifts){ - - int i; - int noerr=1; - - /*output: */ - int* riftsnumpenaltypairs=NULL; - double** riftspenaltypairs=NULL; - int* riftsfill=NULL; - double* riftsfriction=NULL; - - /*intermediary: */ - double fill; - double* riftpenaltypairs=NULL; - #ifdef _SERIAL_ - mxArray* pmxa_rifts=NULL; - #endif - - if(numrifts){ - - /*Allocate arrays: */ - riftsnumpenaltypairs=(int*)xmalloc(numrifts*sizeof(int)); - riftspenaltypairs=(double**)xmalloc(numrifts*sizeof(double*)); - riftsfill=(int*)xmalloc(numrifts*sizeof(int)); - riftsfriction=(double*)xmalloc(numrifts*sizeof(double)); - - #ifdef _SERIAL_ - /*From model handle, recover rifts handle: */ - pmxa_rifts=mxGetField(model_handle,0,"rifts"); - - for(i=0;i +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../shared/shared.h" +#include "../include/typedefs.h" +#include "../include/macros.h" +#include "./objects.h" + + +Riftfront::Riftfront(){ + return; +} + +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_penalty_offset, bool riftfront_penalty_lock, bool riftfront_active,int riftfront_counter,bool riftfront_prestable,bool riftfront_shelf){ + + int i; + + strcpy(type,riftfront_type); + id=riftfront_id; + + for(i=0;iEcho(); + } + printf("\n"); + + printf(" mparid: %i\n",mparid); + if(matpar)matpar->Echo(); + + printf("normal [%g,%g], length %g\n",normal[0],normal[1],normal[2]); + printf("fill: %i friction %g\n",fill,friction); + printf("penalty_offset %g\n",penalty_offset); + printf("penalty_lock %i\n",penalty_lock); + printf("active %i\n",active); + printf("counter %i\n",counter); + printf("prestable %i\n",prestable); + printf("shelf %i\n",shelf); + +} + +void Riftfront::Marshall(char** pmarshalled_dataset){ + + char* marshalled_dataset=NULL; + int enum_type=0; + + /*recover marshalled_dataset: */ + marshalled_dataset=*pmarshalled_dataset; + + /*get enum type of Riftfront: */ + enum_type=RiftfrontEnum(); + + /*marshall enum: */ + memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type); + + /*marshall Riftfront data: */ + memcpy(marshalled_dataset,&type,sizeof(type));marshalled_dataset+=sizeof(type); + memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id); + memcpy(marshalled_dataset,&node_ids,sizeof(node_ids));marshalled_dataset+=sizeof(node_ids); + memcpy(marshalled_dataset,&node_offsets,sizeof(node_offsets));marshalled_dataset+=sizeof(node_offsets); + memcpy(marshalled_dataset,&mparid,sizeof(mparid));marshalled_dataset+=sizeof(mparid); + memcpy(marshalled_dataset,&matpar_offset,sizeof(matpar_offset));marshalled_dataset+=sizeof(matpar_offset); + + memcpy(marshalled_dataset,&h,sizeof(h));marshalled_dataset+=sizeof(h); + memcpy(marshalled_dataset,&b,sizeof(b));marshalled_dataset+=sizeof(b); + memcpy(marshalled_dataset,&s,sizeof(s));marshalled_dataset+=sizeof(s); + + memcpy(marshalled_dataset,&normal,sizeof(normal));marshalled_dataset+=sizeof(normal); + memcpy(marshalled_dataset,&length,sizeof(length));marshalled_dataset+=sizeof(length); + memcpy(marshalled_dataset,&fill,sizeof(fill));marshalled_dataset+=sizeof(fill); + memcpy(marshalled_dataset,&friction,sizeof(friction));marshalled_dataset+=sizeof(friction); + memcpy(marshalled_dataset,&penalty_offset,sizeof(penalty_offset));marshalled_dataset+=sizeof(penalty_offset); + memcpy(marshalled_dataset,&penalty_lock,sizeof(penalty_lock));marshalled_dataset+=sizeof(penalty_lock); + memcpy(marshalled_dataset,&active,sizeof(active));marshalled_dataset+=sizeof(active); + memcpy(marshalled_dataset,&counter,sizeof(counter));marshalled_dataset+=sizeof(counter); + memcpy(marshalled_dataset,&prestable,sizeof(prestable));marshalled_dataset+=sizeof(prestable); + memcpy(marshalled_dataset,&shelf,sizeof(shelf));marshalled_dataset+=sizeof(shelf); + + *pmarshalled_dataset=marshalled_dataset; + return; +} + +int Riftfront::MarshallSize(){ + + return sizeof(type)+ + sizeof(id)+ + sizeof(node_ids)+ + sizeof(node_offsets)+ + sizeof(mparid)+ + sizeof(matpar_offset)+ + sizeof(h)+ + sizeof(b)+ + sizeof(s)+ + sizeof(normal)+ + sizeof(length)+ + sizeof(fill)+ + sizeof(friction)+ + sizeof(penalty_offset)+ + sizeof(penalty_lock)+ + sizeof(active)+ + sizeof(counter)+ + sizeof(prestable)+ + sizeof(shelf)+ + sizeof(int); //sizeof(int) for enum type +} + +char* Riftfront::GetName(void){ + return "riftfront"; +} + + +void Riftfront::Demarshall(char** pmarshalled_dataset){ + + int i; + char* marshalled_dataset=NULL; + + /*recover marshalled_dataset: */ + marshalled_dataset=*pmarshalled_dataset; + + /*this time, no need to get enum type, the pointer directly points to the beginning of the + *object data (thanks to DataSet::Demarshall):*/ + + memcpy(&type,marshalled_dataset,sizeof(type));marshalled_dataset+=sizeof(type); + memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id); + memcpy(&node_ids,marshalled_dataset,sizeof(node_ids));marshalled_dataset+=sizeof(node_ids); + memcpy(&node_offsets,marshalled_dataset,sizeof(node_offsets));marshalled_dataset+=sizeof(node_offsets); + memcpy(&mparid,marshalled_dataset,sizeof(mparid));marshalled_dataset+=sizeof(mparid); + memcpy(&matpar_offset,marshalled_dataset,sizeof(matpar_offset));marshalled_dataset+=sizeof(matpar_offset); + + memcpy(&h,marshalled_dataset,sizeof(h));marshalled_dataset+=sizeof(h); + memcpy(&b,marshalled_dataset,sizeof(b));marshalled_dataset+=sizeof(b); + memcpy(&s,marshalled_dataset,sizeof(s));marshalled_dataset+=sizeof(s); + memcpy(&normal,marshalled_dataset,sizeof(normal));marshalled_dataset+=sizeof(normal); + memcpy(&length,marshalled_dataset,sizeof(length));marshalled_dataset+=sizeof(length); + memcpy(&fill,marshalled_dataset,sizeof(fill));marshalled_dataset+=sizeof(fill); + memcpy(&friction,marshalled_dataset,sizeof(friction));marshalled_dataset+=sizeof(friction); + memcpy(&penalty_offset,marshalled_dataset,sizeof(penalty_offset));marshalled_dataset+=sizeof(penalty_offset); + memcpy(&penalty_lock,marshalled_dataset,sizeof(penalty_lock));marshalled_dataset+=sizeof(penalty_lock); + memcpy(&active,marshalled_dataset,sizeof(active));marshalled_dataset+=sizeof(active); + memcpy(&counter,marshalled_dataset,sizeof(counter));marshalled_dataset+=sizeof(counter); + memcpy(&prestable,marshalled_dataset,sizeof(prestable));marshalled_dataset+=sizeof(prestable); + memcpy(&shelf,marshalled_dataset,sizeof(shelf));marshalled_dataset+=sizeof(shelf); + + for(i=0;iRecover("thickness",&h[0],1,dofs,MAX_RIFTFRONT_GRIDS,(void**)nodes); + inputs->Recover("bed",&b[0],1,dofs,MAX_RIFTFRONT_GRIDS,(void**)nodes); + inputs->Recover("surface",&s[0],1,dofs,MAX_RIFTFRONT_GRIDS,(void**)nodes); + +} + +#undef __FUNCT__ +#define __FUNCT__ "Riftfront::GetDofList" + +void Riftfront::GetDofList(int* doflist,int* pnumberofdofspernode){ + + int i,j; + int doflist_per_node[MAXDOFSPERNODE]; + int numberofdofspernode; + + for(i=0;iGetDofList(&doflist_per_node[0],&numberofdofspernode); + for(j=0;jactive){ + + /*There is contact, we need to constrain the normal velocities (zero penetration), and the + *contact slip friction. */ + + #ifdef _DEBUG_ + printf("Dealing with grid pair (%i,%i)\n",nodes[0]->GetId(),nodes[1]->GetId()); + #endif + + /*Recover input parameters: */ + inputs->Recover("thickness",&h[0],1,dofs,MAX_RIFTFRONT_GRIDS,(void**)nodes); + if (h[0]!=h[1])throw ErrorException(__FUNCT__," different thicknesses not supported for rift fronts"); + thickness=h[0]; + + #ifdef _DEBUG_ + printf("Thickness at grid (%i,%i): %lg\n",nodes[0]->GetId(),nodes[1]->GetID(),thickness); + #endif + + /*From Peter Wriggers book (Computational Contact Mechanics, p191): */ + //First line: + Ke_gg[0][0]+=pow(normal[0],2)*kmax*pow(10,penalty_offset); + Ke_gg[0][1]+=normal[0]*normal[1]*kmax*pow(10,penalty_offset); + Ke_gg[0][2]+=-pow(normal[0],2)*kmax*pow(10,penalty_offset); + Ke_gg[0][3]+=-normal[0]*normal[1]*kmax*pow(10,penalty_offset); + //Second line: + Ke_gg[1][0]+=normal[0]*normal[1]*kmax*pow(10,penalty_offset); + Ke_gg[1][1]+=pow(normal[1],2)*kmax*pow(10,penalty_offset); + Ke_gg[1][2]+=-normal[0]*normal[1]*kmax*pow(10,penalty_offset); + Ke_gg[1][3]+=-pow(normal[1],2)*kmax*pow(10,penalty_offset); + //Third line: + Ke_gg[2][0]+=-pow(normal[0],2)*kmax*pow(10,penalty_offset); + Ke_gg[2][1]+=-normal[0]*normal[1]*kmax*pow(10,penalty_offset); + Ke_gg[2][2]+=pow(normal[0],2)*kmax*pow(10,penalty_offset); + Ke_gg[2][3]+=normal[0]*normal[1]*kmax*pow(10,penalty_offset); + //Fourth line: + Ke_gg[3][0]+=-normal[0]*normal[1]*kmax*pow(10,penalty_offset); + Ke_gg[3][1]+=-pow(normal[1],2)*kmax*pow(10,penalty_offset); + Ke_gg[3][2]+=normal[0]*normal[1]*kmax*pow(10,penalty_offset); + Ke_gg[3][3]+=pow(normal[1],2)*kmax*pow(10,penalty_offset); + + /*Now take care of the friction: of type sigma=frictiontangent_velocity2-tangent_velocity1)*/ + + //First line: + Ke_gg[0][0]+=pow(normal[1],2)*thickness*length*friction; + Ke_gg[0][1]+=-normal[0]*normal[1]*thickness*length*friction; + Ke_gg[0][2]+=-pow(normal[1],2)*thickness*length*friction; + Ke_gg[0][3]+=normal[0]*normal[1]*thickness*length*friction; + //Second line: + Ke_gg[1][0]+=-normal[0]*normal[1]*thickness*length*friction; + Ke_gg[1][1]+=pow(normal[0],2)*thickness*length*friction; + Ke_gg[1][2]+=normal[0]*normal[1]*thickness*length*friction; + Ke_gg[1][3]+=-pow(normal[0],2)*thickness*length*friction; + //Third line: + Ke_gg[2][0]+=-pow(normal[1],2)*thickness*length*friction; + Ke_gg[2][1]+=normal[0]*normal[1]*thickness*length*friction; + Ke_gg[2][2]+=pow(normal[1],2)*thickness*length*friction; + Ke_gg[2][3]+=-normal[0]*normal[1]*thickness*length*friction; + //Fourth line: + Ke_gg[3][0]+=normal[0]*normal[1]*thickness*length*friction; + Ke_gg[3][1]+=-pow(normal[0],2)*thickness*length*friction; + Ke_gg[3][2]+=-normal[0]*normal[1]*thickness*length*friction; + Ke_gg[3][3]+=pow(normal[0],2)*thickness*length*friction; + + /*Add Ke_gg to global matrix Kgg: */ + MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES); + } + else{ + /*the grids on both sides of the rift do not penetrate. PenaltyCreatePVector will + *take care of adding point loads to simulate pressure on the rift flanks. But as far as stiffness, + there is none (0 stiffness implies decoupling of the flank rifts, which is exactly what we want): */ + } + +} + +#undef __FUNCT__ +#define __FUNCT__ "Riftfront::PenaltyCreatePVector" +void Riftfront::PenaltyCreatePVector(Vec pg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type){ + + int i,j; + const int numgrids=MAX_RIFTFRONT_GRIDS; + int dofs[1]={0}; + double pe_g[4]; + const int numdof=2*numgrids; + int doflist[numdof]; + int numberofdofspernode; + ParameterInputs* inputs=NULL; + double rho_ice; + double rho_water; + double gravity; + double thickness; + double bed; + double pressure; + + /*Some pointer intialization: */ + inputs=(ParameterInputs*)vinputs; + + /* Get node coordinates and dof list: */ + GetDofList(&doflist[0],&numberofdofspernode); + + /* Set pe_g to 0: */ + for(i=0;iactive){ + /*Ok, this rift is opening. We should put loads on both sides of the rift flanks. Because we are dealing with contact mechanics, + * and we want to avoid zigzagging of the loads, we want lump the loads onto grids, not onto surfaces between grids.:*/ + + #ifdef _DEBUG_ + _printf_("Grids (%i,%i) are free of constraints\n",nodes[0]->GetId(),nodes[1]->GetID()); + #endif + + /*Ok, to compute the pressure, we are going to need material properties, thickness and bed for the two grids. We assume those properties to + * be the same across the rift.: */ + + rho_ice=matpar->GetRhoIce(); + rho_water=matpar->GetRhoWater(); + gravity=matpar->GetG(); + + /*get thickness: */ + inputs->Recover("thickness",&h[0],1,dofs,MAX_RIFTFRONT_GRIDS,(void**)nodes); + if (h[0]!=h[1])throw ErrorException(__FUNCT__," different thicknesses not supported for rift fronts"); + thickness=h[0]; + + inputs->Recover("bed",&b[0],1,dofs,MAX_RIFTFRONT_GRIDS,(void**)nodes); + if (b[0]!=b[1])throw ErrorException(__FUNCT__," different beds not supported for rift fronts"); + bed=b[0]; + + /*Ok, now compute the pressure (in norm) that is being applied to the flanks, depending on the type of fill: */ + if(fill==WATERFILL){ + if(shelf){ + /*We are on an ice shelf, hydrostatic equilibrium is used to determine the pressure for water fill: */ + pressure=rho_ice*gravity*pow(thickness,(double)2)/(double)2 - rho_water*gravity*pow(bed,(double)2)/(double)2; + } + else{ + //We are on an icesheet, we assume the water column fills the entire front: */ + pressure=rho_ice*gravity*pow(thickness,(double)2)/(double)2 - rho_water*gravity*pow(thickness,(double)2)/(double)2; + } + } + else if(fill==AIRFILL){ + pressure=rho_ice*gravity*pow(thickness,(double)2)/(double)2; //icefront on an ice sheet, pressure imbalance ice vs air. + } + else if(fill==ICEFILL){ //icefront finding itself against another icefront (pressure imbalance is fully compensated, ice vs ice) + pressure=0; + } + else{ + throw ErrorException(__FUNCT__,exprintf("%s%s%i%s",__FUNCT__," error message: fill type ",fill," not supported yet.")); + } + + /*Ok, add contribution to first grid, along the normal i==0: */ + for (j=0;j<2;j++){ + pe_g[j]+=pressure*normal[j]*length; + } + + /*Add contribution to second grid, along the opposite normal: i==1 */ + for (j=0;j<2;j++){ + pe_g[2+j]+= -pressure*normal[j]*length; + } + /*Add pe_g to global vector pg; */ + VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES); + + } + else{ + /*The penalty is active. No loads implied here.*/ + } +} + +Object* Riftfront::copy() { + return new Riftfront(*this); +} + +#undef __FUNCT__ +#define __FUNCT__ "Riftfront::CreateKMatrix" +void Riftfront::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){ + /*do nothing: */ +} + +#undef __FUNCT__ +#define __FUNCT__ "Riftfront::CreatePVector" +void Riftfront::CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type){ + /*do nothing: */ +} + +bool Riftfront::PreStable(){ + return prestable; +} + +void Riftfront::SetPreStable(){ + prestable=1; +} + +#undef __FUNCT__ +#define __FUNCT__ "Riftfront::PreConstrain" +int Riftfront::PreConstrain(int* punstable, void* vinputs, int analysis_type){ + + const int numgrids=2; + int dofs[2]={0,1}; + double vxvy_list[2][2]; //velocities for all grids + double penetration; + int unstable; + ParameterInputs* inputs=NULL; + int found; + + inputs=(ParameterInputs*)vinputs; + + /*First recover velocity: */ + found=inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes); + if(!found)throw ErrorException(__FUNCT__," could not find velocity in inputs!"); + + /*Grid 1 faces grid2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */ + penetration=(vxvy_list[1][0]-vxvy_list[0][0])*normal[0]+(vxvy_list[1][1]-vxvy_list[0][1])*normal[1]; + + /*Ok, we are preconstraining here. Ie, anything that penetrates is constrained until stability of the entire set + * of constraints is reached.: */ + if(penetration<0){ + if (!this->active){ + /*This is the first time penetration happens: */ + this->active=1; + unstable=1; + } + else{ + /*This constraint was already active: */ + this->active=1; + unstable=0; + } + } + else{ + /*No penetration happening. : */ + if (!this->active){ + /*This penalty was not active, and no penetration happening. Do nonthing: */ + this->active=0; + unstable=0; + } + else{ + /*Ok, this penalty wants to get released. But not now, this is preconstraint, not constraint: */ + this->active=1; + unstable=0; + } + } + + /*assign output pointer: */ + *punstable=unstable; +} + +#undef __FUNCT__ +#define __FUNCT__ "Riftfront::Constrain" +int Riftfront::Constrain(int* punstable, void* vinputs, int analysis_type){ + + const int numgrids=2; + int dofs[2]={0,1}; + double vxvy_list[2][2]; //velocities for all grids + double max_penetration; + double penetration; + int activate; + int found; + int unstable; + + ParameterInputs* inputs=NULL; + + inputs=(ParameterInputs*)vinputs; + + + /*First recover parameter inputs: */ + found=inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes); + if(!found)throw ErrorException(__FUNCT__," could not find velocity in inputs!"); + + found=inputs->Recover("max_penetration",&max_penetration); + if(!found)throw ErrorException(__FUNCT__," could not find max_penetration in inputs!"); + + /*Grid 1 faces grid2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */ + penetration=(vxvy_list[1][0]-vxvy_list[0][0])*normal[0]+(vxvy_list[1][1]-vxvy_list[0][1])*normal[1]; + + /*Activate or deactivate penalties: */ + if(penetration<0){ + /*There is penetration, we need to active the penalty so this penetration will be NULL: */ + activate=1; + } + else{ + /*Ok, there is no penetration for these two grids of the rift. We want to deactivate the penalty. If we do + * it all at once, then we will zigzag forever. Only deactivate once at a time. Deactivate the one that most + * wants to, ie the one with the max penetration: */ + if (penetration==max_penetration){ + activate=0; + } + else{ + /*Ok, we are dealing with the rest of the penalties that want to be freed. But may be in this lot there + * are penalties that were already free? Keep them as such. For the rest, do not release them yet: */ + if (this->active==0){ + activate=0; + } + else{ + activate=1; + } + } + } + + /*Here, we try to avoid zigzaging. When a penalty activates and deactivates for more than penalty_lock times, + * we lock it: */ + if(this->counter>this->penalty_lock){ + /*This penalty has zig zagged too long, fix it to smooth results: */ + activate=1; this->active=activate; + printf(" locking riftfront %i\n",this->id); + } + + //Figure out stability of this penalty + if(this->active==activate){ + unstable=0; + } + else{ + unstable=1; + this->counter++; + } + + //Set penalty flag + this->active=activate; + + /*assign output pointer: */ + *punstable=unstable; +} + +#undef __FUNCT__ +#define __FUNCT__ "Riftfront::Penetration" +int Riftfront::Penetration(double* ppenetration, void* vinputs, int analysis_type){ + + const int numgrids=2; + int dofs[2]={0,1}; + double vxvy_list[2][2]; //velocities for all grids + double max_penetration; + double penetration; + int found; + + ParameterInputs* inputs=NULL; + + inputs=(ParameterInputs*)vinputs; + + + found=inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes); + if(!found)throw ErrorException(__FUNCT__," could not find velocity in inputs!"); + + /*Grid 1 faces grid2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */ + penetration=(vxvy_list[1][0]-vxvy_list[0][0])*normal[0]+(vxvy_list[1][1]-vxvy_list[0][1])*normal[1]; + + /*Now, we return penetration only if we are active!: */ + if(this->active==0)penetration=0; + + /*assign output pointer: */ + *ppenetration=penetration; + +} + +#undef __FUNCT__ +#define __FUNCT__ "Riftfront::PotentialUnstableConstraint" +int Riftfront::PotentialUnstableConstraint(int* punstable, void* vinputs, int analysis_type){ + + + const int numgrids=2; + int dofs[2]={0,1}; + double vxvy_list[2][2]; //velocities for all grids + double max_penetration; + double penetration; + int activate; + int unstable; + int found; + + ParameterInputs* inputs=NULL; + + inputs=(ParameterInputs*)vinputs; + + found=inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes); + if(!found)throw ErrorException(__FUNCT__," could not find velocity in inputs!"); + + /*Grid 1 faces grid2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */ + penetration=(vxvy_list[1][0]-vxvy_list[0][0])*normal[0]+(vxvy_list[1][1]-vxvy_list[0][1])*normal[1]; + + /*Ok, we are looking for positive penetration in an active constraint: */ + if(this->active){ + if (penetration>=0){ + unstable=1; + } + else{ + unstable=0; + } + } + else{ + unstable=0; + } + + /*assign output pointer: */ + *punstable=unstable; +} Index: /issm/trunk/src/c/objects/Riftfront.h =================================================================== --- /issm/trunk/src/c/objects/Riftfront.h (revision 1628) +++ /issm/trunk/src/c/objects/Riftfront.h (revision 1628) @@ -0,0 +1,90 @@ +/*!\file Riftfront.h + * \brief: header file for riftfront object + */ + +#ifndef _RIFTFRONT_H_ +#define _RIFTFRONT_H_ + +#include "./Load.h" +#include "./Matpar.h" +#include "./Element.h" +#include "./Node.h" +#include "./ParameterInputs.h" + +#define MAX_RIFTFRONT_GRIDS 2 //max number of grids on a rift flank, only 2 because 2d for now. +#define RIFTFRONTSTRING 20 //max string length + +/*Types of rift filling: */ +#define WATERFILL 1 +#define AIRFILL 2 +#define ICEFILL 3 + +class Element; +class Riftfront: public Load { + + private: + char type[RIFTFRONTSTRING]; + int id; + + /*nodes: */ + int node_ids[MAX_RIFTFRONT_GRIDS]; //node ids + Node* nodes[MAX_RIFTFRONT_GRIDS]; //node pointers + int node_offsets[MAX_RIFTFRONT_GRIDS]; //node offsets in nodes dataset + + /*material: */ + int mparid; + Matpar* matpar; + int matpar_offset; + + /*properties: */ + double h[MAX_RIFTFRONT_GRIDS]; //thickness + double b[MAX_RIFTFRONT_GRIDS]; //bed + double s[MAX_RIFTFRONT_GRIDS]; //surface + + double normal[2]; + double length; + int fill; + double friction; + bool shelf; + + double penalty_offset; + bool penalty_lock; + + /*computational: */ + bool active; + int counter; + bool prestable; + + + public: + + Riftfront(); + 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 penalty_offset, bool penalty_lock,bool active,int counter,bool prestable,bool shelf); + ~Riftfront(); + + void Echo(); + void DeepEcho(); + void Marshall(char** pmarshalled_dataset); + int MarshallSize(); + char* GetName(); + void Demarshall(char** pmarshalled_dataset); + int Enum(); + int GetId(); + int MyRank(); + void Configure(void* elements,void* nodes,void* materials); + void UpdateFromInputs(void* inputs); + void GetDofList(int* doflist,int* pnumberofdofs); + void CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type); + void CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type); + 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); + bool PreStable(); + void SetPreStable(); + int PreConstrain(int* punstable, void* inputs, int analysis_type); + int Constrain(int* punstable, void* inputs, int analysis_type); + int Penetration(double* ppenetration, void* inputs, int analysis_type); + int PotentialUnstableConstraint(int* punstable, void* inputs, int analysis_type); + Object* copy(); +}; + +#endif /* _RIFTFRONT_H_ */ Index: /issm/trunk/src/c/objects/objects.h =================================================================== --- /issm/trunk/src/c/objects/objects.h (revision 1627) +++ /issm/trunk/src/c/objects/objects.h (revision 1628) @@ -22,4 +22,5 @@ #include "./Rgb.h" #include "./Icefront.h" +#include "./Riftfront.h" #include "./Penpair.h" #include "./Pengrid.h"