Changeset 27489

Timestamp:

01/01/23 07:32:41 (2 years ago)

Author:

Mathieu Morlighem

Message:

CHG: fone with linearize=2 (elementwise friction)

File:

• issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (modified) (7 diffs)

issm/trunk-jpl/src/c/classes/Loads/Friction.cpp

@@ -34 +34 @@
          * There are exceptions, e.g. HO, which needs the user to specify the dimension used in Friciton.*/
 
+        /*Intermediaries*/
+        int linearization_type;
+
         this->element=element_in;
         this->linearize  = 0;
@@ -65 +68 @@
 
         if(this->law==1){
-                element_in->FindParam(&this->linearize,FrictionLinearizeEnum);
-                if(this->linearize){
-                        this->linearize = 0; /*Change to make sure we do the calculation once*/
+                element_in->FindParam(&linearization_type,FrictionLinearizeEnum);
+                if(linearization_type==0){
+                        /*Don't do anything*/
+                }
+                else if(linearization_type==1){
                         int numvertices = this->element->GetNumberOfVertices();
                         this->alpha2_list            = xNew<IssmDouble>(numvertices);
@@ -75 +80 @@
                                 gauss->GaussVertex(iv);
                                 this->GetAlpha2(&this->alpha2_list[iv], gauss);
-                                IssmDouble temp = 0.;
-                                this->GetAlpha2(&temp, gauss);
                                 this->GetAlphaComplement(&this->alpha2_complement_list[iv], gauss);
                         }
-                        this->linearize = 1; /*Change back, we are now all set!*/
+                        this->linearize = linearization_type; /*Change back, we are now all set!*/
+                        delete gauss;
+                }
+                else if(linearization_type==2){
+                        this->alpha2_list            = xNew<IssmDouble>(1);
+                        this->alpha2_complement_list = xNew<IssmDouble>(1);
+                        Gauss* gauss=element->NewGauss(1); gauss->GaussPoint(0);
+                        this->GetAlpha2(&this->alpha2_list[0], gauss);
+                        this->GetAlphaComplement(&this->alpha2_complement_list[0], gauss);
+                        this->linearize = linearization_type; /*Change back, we are now all set!*/
+                        delete gauss;
+                }
+                else{
+                        _error_("not supported yet");
                 }
         }
@@ -94 +110 @@
 Friction::~Friction(){/*{{{*/
         if(this->linearize){
-                printf("linearize = %i\n",this->linearize);
                 xDelete<IssmDouble>(this->alpha2_list);
                 xDelete<IssmDouble>(this->alpha2_complement_list);
@@ -111 +126 @@
         switch(this->law){
                 case 1:
-                        GetAlphaViscousComplement(palpha_complement,gauss);
+                        if(this->linearize==0){
+                                GetAlphaViscousComplement(palpha_complement,gauss);
+                        }
+                        else if(this->linearize==1){
+                                this->element->ValueP1OnGauss(palpha_complement, this->alpha2_complement_list, gauss);
+                        }
+                        else if(this->linearize==2){
+                                *palpha_complement = this->alpha2_complement_list[0];
+                        }
+                        else{
+                                _error_("not supported yet");
+                        }
+                        break;
                         break;
                 case 2:
@@ -333 +360 @@
         switch(this->law){
                 case 1:
-                        if(this->linearize){
+                        if(this->linearize==0){
+                                GetAlpha2Viscous(palpha2,gauss);
+                        }
+                        else if(this->linearize==1){
                                 this->element->ValueP1OnGauss(palpha2, this->alpha2_list, gauss);
                         }
+                        else if(this->linearize==2){
+                                *palpha2 = this->alpha2_list[0];
+                        }
                         else{
-                                GetAlpha2Viscous(palpha2,gauss);
+                                _error_("not supported yet");
                         }
                         break;
@@ -1288 +1321 @@
         switch(frictionlaw){
                 case 1:
-                        //parameters->AddObject(iomodel->CopyConstantObject("md.friction.linearize",FrictionLinearizeEnum));
+                        parameters->AddObject(iomodel->CopyConstantObject("md.friction.linearize",FrictionLinearizeEnum));
                         parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
                         parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum));