Changeset 4978

Timestamp:

08/04/10 12:09:21 (15 years ago)

Author:

Mathieu Morlighem

Message:

Enable CM RheologyB2d in 3d, still some work to do

Location:

issm/trunk/src

Files:

• c/EnumDefinitions/EnumAsString.cpp (modified) (1 diff)
• c/EnumDefinitions/EnumDefinitions.h (modified) (1 diff)
• c/EnumDefinitions/StringAsEnum.cpp (modified) (1 diff)
• c/modules/ModelProcessorx/Control/CreateParametersControl.cpp (modified) (1 diff)
• c/objects/Elements/Penta.cpp (modified) (9 diffs)
• c/objects/Elements/Penta.h (modified) (1 diff)
• c/objects/Materials/Matice.cpp (modified) (1 diff)
• m/enum/ExtrudeParamEnum.m (deleted)

issm/trunk/src/c/EnumDefinitions/EnumAsString.cpp

@@ -268 +268 @@
                 case EpsRelEnum : return "EpsRel";
                 case EpsResEnum : return "EpsRes";
-                case ExtrudeParamEnum : return "ExtrudeParam";
                 case HeatCapacityEnum : return "HeatCapacity";
                 case IsHutterEnum : return "IsHutter";

issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h

@@ -305 +305 @@
         EpsRelEnum,
         EpsResEnum,
-        ExtrudeParamEnum,
         HeatCapacityEnum,
         IsHutterEnum,

issm/trunk/src/c/EnumDefinitions/StringAsEnum.cpp

@@ -266 +266 @@
         else if (strcmp(name,"EpsRel")==0) return EpsRelEnum;
         else if (strcmp(name,"EpsRes")==0) return EpsResEnum;
-        else if (strcmp(name,"ExtrudeParam")==0) return ExtrudeParamEnum;
         else if (strcmp(name,"HeatCapacity")==0) return HeatCapacityEnum;
         else if (strcmp(name,"IsHutter")==0) return IsHutterEnum;

issm/trunk/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp

@@ -26 +26 @@
 
                 if (strcmp(iomodel->control_type,"drag_coefficient")==0){
-                        parameters->AddObject(new BoolParam(ExtrudeParamEnum,false));
                         parameters->AddObject(new IntParam(ControlTypeEnum,DragCoefficientEnum));
                 }
                 else  if (strcmp(iomodel->control_type,"rheology_B")==0){
-                        parameters->AddObject(new BoolParam(ExtrudeParamEnum,true));
                         parameters->AddObject(new IntParam(ControlTypeEnum,RheologyB2dEnum));
                 }

issm/trunk/src/c/objects/Elements/Penta.cpp

@@ -697 +697 @@
                 /*This element should be collapsed into a tria element at its base. Create this tria element, 
                  * and compute CostFunction*/
+
+                /*Depth Average B*/
+                this->InputDepthAverageAtBase(RheologyBEnum,RheologyB2dEnum,MaterialsEnum);
+
                 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
                 J=tria->CostFunction(process_units);
                 delete tria;
+
+                /*delete B average*/
+                this->matice->inputs->DeleteInput(RheologyB2dEnum);
+
                 return J;
         }
         else{
+
+                /*Depth Average B*/
+                this->InputDepthAverageAtBase(RheologyBEnum,RheologyB2dEnum,MaterialsEnum);
 
                 tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
                 J=tria->CostFunction(process_units);
                 delete tria;
+
+                /*delete B average*/
+                this->matice->inputs->DeleteInput(RheologyB2dEnum);
+
                 return J;
         }
@@ -994 +1009 @@
         double cm_min,cm_max;
         int    control_type;
+        bool   onbed;
 
         /*retrieve some parameters: */
@@ -1000 +1016 @@
         this->parameters->FindParam(&control_type,ControlTypeEnum);
 
-
         /*Rheology*/
         if(control_type==RheologyB2dEnum){
+
+                /*The update is done by the element on bed only*/
+                inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+                if (!onbed) return;
 
                 /*First, get revert to previous parameter value (kept in ControlParameter input)*/
@@ -1018 +1037 @@
                 input->Constrain(cm_min,cm_max);
 
+                /*OK, now we can propagate the new parameter up*/
+                this->InputExtrude(RheologyBEnum,MaterialsEnum);
+
                 /*Finally: save input if requested*/
-                if (save_parameter) matice->inputs->DuplicateInput(RheologyBEnum,ControlParameterEnum);
+                if (save_parameter){
+                        matice->inputs->DuplicateInput(RheologyBEnum,ControlParameterEnum);
+                        this->InputExtrude(ControlParameterEnum,MaterialsEnum);
+                }
 
         }
@@ -1174 +1199 @@
         /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */
         if (enum_type==RheologyB2dEnum){
-                input=this->matice->inputs->GetInput(enum_type);
+                input=this->matice->inputs->GetInput(RheologyBEnum);
         }
         else{
@@ -4671 +4696 @@
                 /*This element should be collapsed into a tria element at its base. Create this tria element, 
                  * and compute gardj*/
+
+                /*Depth Average B*/
+                this->InputDepthAverageAtBase(RheologyBEnum,RheologyB2dEnum,MaterialsEnum);
+
                 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
                 tria->GradjB(gradient);
                 delete tria;
+
+                /*delete B average*/
+                this->matice->inputs->DeleteInput(RheologyB2dEnum);
         }
         else{
+                /*Depth Average B*/
+                this->InputDepthAverageAtBase(RheologyBEnum,RheologyB2dEnum,MaterialsEnum);
+
                 /*B is a 2d field, use MacAyeal(2d) gradient even if it is Stokes or Pattyn*/
                 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
                 tria->GradjB(gradient);
                 delete tria;
-        }
-
+
+                /*delete B average*/
+                this->matice->inputs->DeleteInput(RheologyB2dEnum);
+        }
 
 }
@@ -4736 +4773 @@
 /*}}}*/
 /*FUNCTION Penta::InputExtrude {{{1*/
-void  Penta::InputExtrude(int enum_type,bool only_if_collapsed){
-
-        bool   onbed,collapse=false;
+void  Penta::InputExtrude(int enum_type,int object_type){
+
+        bool   onbed;
         Penta *penta          = NULL;
         Input *original_input = NULL;
 
         /*recover parameters: */
-        if (only_if_collapsed) inputs->GetParameterValue(&collapse,CollapseEnum);
         inputs->GetParameterValue(&onbed,ElementOnBedEnum);
-
-        /*First: if only_if_collapsed, check wether this penta is collapsed*/
-        if (only_if_collapsed && !collapse) return;
 
         /*Are we on the base, not on the surface?:*/
@@ -4758 +4791 @@
 
                 /*Step1: Extrude the original input: */
-                original_input=(Input*)this->inputs->GetInput(enum_type);
+                if (object_type==ElementsEnum)
+                 original_input=(Input*)this->inputs->GetInput(enum_type);
+                else if (object_type==MaterialsEnum)
+                 original_input=(Input*)matice->inputs->GetInput(enum_type);
+                else
+                 ISSMERROR("object of type %s not supported yet",EnumAsString(object_type));
                 if(!original_input) ISSMERROR("%s%s"," could not find input with enum:",EnumAsString(enum_type));
                 original_input->Extrude();
@@ -4776 +4814 @@
                         Input* copy=NULL;
                         copy=(Input*)original_input->copy();
-                        penta->inputs->AddInput((Input*)copy);
+                        if (object_type==ElementsEnum)
+                         penta->inputs->AddInput((Input*)copy);
+                        else if (object_type==MaterialsEnum)
+                         penta->matice->inputs->AddInput((Input*)copy);
+                        else
+                         ISSMERROR("object of type %s not supported yet",EnumAsString(object_type));
 
                         /*Stop if we have reached the surface*/

issm/trunk/src/c/objects/Elements/Penta.h

@@ -156 +156 @@
                 void    GetStrainRate3d(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input, Input* vz_input);
                 Penta*  GetUpperElement(void);
-                void      InputExtrude(int enum_type,bool only_if_collapsed);
+                void      InputExtrude(int enum_type,int object_type);
                 void    InputUpdateFromSolutionAdjointHoriz( double* solutiong);
                 void    InputUpdateFromSolutionAdjointStokes( double* solutiong);

issm/trunk/src/c/objects/Materials/Matice.cpp

@@ -547 +547 @@
 
         /*Call inputs method*/
-        printf("%s\n",EnumAsString(original_enum));
         if (IsInput(original_enum)) inputs->DuplicateInput(original_enum,new_enum);