[issm-svn] r19127 - issm/trunk/src/c/analyses
morlighe at issm.ess.uci.edu
morlighe at issm.ess.uci.edu
Thu Feb 19 16:23:50 PST 2015
Author: morlighe
Date: 2015-02-19 16:23:50 -0800 (Thu, 19 Feb 2015)
New Revision: 19127
Modified:
issm/trunk/src/c/analyses/LevelsetAnalysis.cpp
Log:
BUG: LevelsetAnalysis.cpp was not updated properly
Modified: issm/trunk/src/c/analyses/LevelsetAnalysis.cpp
===================================================================
--- issm/trunk/src/c/analyses/LevelsetAnalysis.cpp 2015-02-19 23:48:11 UTC (rev 19126)
+++ issm/trunk/src/c/analyses/LevelsetAnalysis.cpp 2015-02-20 00:23:50 UTC (rev 19127)
@@ -10,19 +10,28 @@
#include "../modules/modules.h"
#include "../solutionsequences/solutionsequences.h"
-int LevelsetAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
- return 1;
+void LevelsetAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
+ return;
}
/*}}}*/
-void LevelsetAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
+void LevelsetAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
return;
+}/*}}}*/
+void LevelsetAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
+ int finiteelement=P1Enum;
+ if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
+ ::CreateNodes(nodes,iomodel,LevelsetAnalysisEnum,finiteelement);
+ iomodel->DeleteData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
}
/*}}}*/
+int LevelsetAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
+ return 1;
+}
+/*}}}*/
void LevelsetAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
- int finiteelement;
/*Finite element type*/
- finiteelement = P1Enum;
+ int finiteelement = P1Enum;
/*Update elements: */
int counter=0;
@@ -37,45 +46,71 @@
iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
iomodel->FetchDataToInput(elements,VxEnum);
iomodel->FetchDataToInput(elements,VyEnum);
- iomodel->FetchDataToInput(elements,MasstransportCalvingrateEnum);
+
+ /*Get calving parameters*/
+ bool iscalving;
+ int calvinglaw;
+ iomodel->Constant(&iscalving,TransientIscalvingEnum);
+ if(iscalving){
+ iomodel->Constant(&calvinglaw,CalvingLawEnum);
+ iomodel->Constant(&iscalving,TransientIscalvingEnum);
+ switch(calvinglaw){
+ case DefaultCalvingEnum:
+ iomodel->FetchDataToInput(elements,CalvingCalvingrateEnum);
+ iomodel->FetchDataToInput(elements,CalvingMeltingrateEnum);
+ break;
+ case CalvingLevermannEnum:
+ iomodel->FetchDataToInput(elements,CalvinglevermannCoeffEnum);
+ iomodel->FetchDataToInput(elements,CalvinglevermannMeltingrateEnum);
+ break;
+ case CalvingPiEnum:
+ iomodel->FetchDataToInput(elements,CalvingpiCoeffEnum);
+ iomodel->FetchDataToInput(elements,CalvingpiMeltingrateEnum);
+ break;
+ case CalvingDevEnum:
+ iomodel->FetchDataToInput(elements,CalvingpiCoeffEnum);
+ iomodel->FetchDataToInput(elements,CalvingMeltingrateEnum);
+ break;
+ default:
+ _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
+ }
+ }
}
/*}}}*/
-void LevelsetAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
- int finiteelement=P1Enum;
- if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
- ::CreateNodes(nodes,iomodel,LevelsetAnalysisEnum,finiteelement);
- iomodel->DeleteData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
-}
-/*}}}*/
-void LevelsetAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
+void LevelsetAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
+ parameters->AddObject(iomodel->CopyConstantObject(LevelsetStabilizationEnum));
return;
}
/*}}}*/
-void LevelsetAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
- return;
-}/*}}}*/
/*Finite element Analysis*/
-void LevelsetAnalysis::Core(FemModel* femmodel){/*{{{*/
+void LevelsetAnalysis::Core(FemModel* femmodel){/*{{{*/
#if !defined(_DEVELOPMENT_)
_error_("Not implemented yet");
#endif
/*parameters: */
+ int stabilization;
bool save_results;
femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
+ femmodel->parameters->FindParam(&stabilization,LevelsetStabilizationEnum);
/*activate formulation: */
femmodel->SetCurrentConfiguration(LevelsetAnalysisEnum);
if(VerboseSolution()) _printf0_("call computational core:\n");
- solutionsequence_linear(femmodel);
+ if(stabilization==4){
+ solutionsequence_fct(femmodel);
+ }
+ else{
+ solutionsequence_linear(femmodel);
+ }
if(save_results){
if(VerboseSolution()) _printf0_(" saving results\n");
- int outputs = MaskIceLevelsetEnum;
- femmodel->RequestedOutputsx(&femmodel->results,&outputs,1);
+ int outputs[2] = {MaskIceLevelsetEnum, CalvingCalvingrateEnum};
+ femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);
}
}/*}}}*/
ElementVector* LevelsetAnalysis::CreateDVector(Element* element){/*{{{*/
@@ -96,31 +131,7 @@
if(!element->IsOnBase()) return NULL;
_error_("not implemented yet");
}/*}}}*/
-void LevelsetAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
- _error_("not implemented yet");
-}/*}}}*/
-void LevelsetAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
- _error_("Not implemented yet");
-}/*}}}*/
-void LevelsetAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
-
- int domaintype;
- element->FindParam(&domaintype,DomainTypeEnum);
- switch(domaintype){
- case Domain2DhorizontalEnum:
- element->InputUpdateFromSolutionOneDof(solution,MaskIceLevelsetEnum);
- break;
- case Domain3DEnum:
- element->InputUpdateFromSolutionOneDofCollapsed(solution,MaskIceLevelsetEnum);
- break;
- default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
- }
-}/*}}}*/
-void LevelsetAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
- /*Default, do nothing*/
- return;
-}/*}}}*/
-void LevelsetAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void LevelsetAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
/*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
* For node i, Bi can be expressed in the actual coordinate system
* by:
@@ -147,7 +158,7 @@
/*Clean-up*/
xDelete<IssmDouble>(basis);
}/*}}}*/
-void LevelsetAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void LevelsetAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
/*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
* For node i, Bi' can be expressed in the actual coordinate system
* by:
@@ -175,8 +186,52 @@
xDelete<IssmDouble>(dbasis);
}/*}}}*/
-void LevelsetAnalysis::SetDistanceOnIntersectedElements(FemModel* femmodel){/*{{{*/
+IssmDouble LevelsetAnalysis::GetDistanceToStraight(IssmDouble* q, IssmDouble* s0, IssmDouble* s1){/*{{{*/
+ // returns distance d of point q to straight going through points s0, s1
+ // d=|a x b|/|b|
+ // with a=q-s0, b=s1-s0
+
+ /* Intermediaries */
+ const int dim=2;
+ int i;
+ IssmDouble a[dim], b[dim];
+ IssmDouble norm_b;
+ for(i=0;i<dim;i++){
+ a[i]=q[i]-s0[i];
+ b[i]=s1[i]-s0[i];
+ }
+
+ norm_b=0.;
+ for(i=0;i<dim;i++)
+ norm_b+=b[i]*b[i];
+ norm_b=sqrt(norm_b);
+ _assert_(norm_b>0.);
+
+ return fabs(a[0]*b[1]-a[1]*b[0])/norm_b;
+}/*}}}*/
+void LevelsetAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
+ _error_("not implemented yet");
+}/*}}}*/
+void LevelsetAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
+ _error_("Not implemented yet");
+}/*}}}*/
+void LevelsetAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
+
+ int domaintype;
+ element->FindParam(&domaintype,DomainTypeEnum);
+ switch(domaintype){
+ case Domain2DhorizontalEnum:
+ element->InputUpdateFromSolutionOneDof(solution,MaskIceLevelsetEnum);
+ break;
+ case Domain3DEnum:
+ element->InputUpdateFromSolutionOneDofCollapsed(solution,MaskIceLevelsetEnum);
+ break;
+ default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
+ }
+}/*}}}*/
+void LevelsetAnalysis::SetDistanceOnIntersectedElements(FemModel* femmodel){/*{{{*/
+
/* Intermediaries */
int i,k;
IssmDouble dmaxp=0.,dmaxm=0,val=0.;
@@ -186,11 +241,11 @@
Element* element;
Vector<IssmDouble>* vec_dist_zerolevelset = NULL;
- GetVectorFromInputsx(&vec_dist_zerolevelset, femmodel, MaskIceLevelsetEnum, VertexEnum);
+ GetVectorFromInputsx(&vec_dist_zerolevelset, femmodel, MaskIceLevelsetEnum, VertexPIdEnum);
/* set NaN on elements intersected by zero levelset */
for(i=0;i<femmodel->elements->Size();i++){
- element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+ element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
if(element->IsZeroLevelset(MaskIceLevelsetEnum))
for(k=0;k<element->GetNumberOfVertices();k++)
vec_dist_zerolevelset->SetValue(element->vertices[k]->Sid(),NAN,INS_VAL);
@@ -198,7 +253,7 @@
/* set distance on elements intersected by zero levelset */
for(i=0;i<femmodel->elements->Size();i++){
- element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+ element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
if(element->IsZeroLevelset(MaskIceLevelsetEnum)){
SetDistanceToZeroLevelsetElement(vec_dist_zerolevelset, element);
@@ -231,7 +286,7 @@
delete vec_dist_zerolevelset;
delete dist_zerolevelset;
}/*}}}*/
-void LevelsetAnalysis::SetDistanceToZeroLevelsetElement(Vector<IssmDouble>* vec_signed_dist, Element* element){/*{{{*/
+void LevelsetAnalysis::SetDistanceToZeroLevelsetElement(Vector<IssmDouble>* vec_signed_dist, Element* element){/*{{{*/
if(!element->IsZeroLevelset(MaskIceLevelsetEnum))
return;
@@ -284,27 +339,7 @@
xDelete<IssmDouble>(xyz_list);
xDelete<IssmDouble>(xyz_list_zero);
}/*}}}*/
-IssmDouble LevelsetAnalysis::GetDistanceToStraight(IssmDouble* q, IssmDouble* s0, IssmDouble* s1){/*{{{*/
- // returns distance d of point q to straight going through points s0, s1
- // d=|a x b|/|b|
- // with a=q-s0, b=s1-s0
-
- /* Intermediaries */
- const int dim=2;
- int i;
- IssmDouble a[dim], b[dim];
- IssmDouble norm_b;
-
- for(i=0;i<dim;i++){
- a[i]=q[i]-s0[i];
- b[i]=s1[i]-s0[i];
- }
-
- norm_b=0.;
- for(i=0;i<dim;i++)
- norm_b+=b[i]*b[i];
- norm_b=sqrt(norm_b);
- _assert_(norm_b>0.);
-
- return fabs(a[0]*b[1]-a[1]*b[0])/norm_b;
+void LevelsetAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
+ /*Default, do nothing*/
+ return;
}/*}}}*/
More information about the issm-svn
mailing list