Changeset 26362

Timestamp:

07/26/21 14:15:53 (4 years ago)

Author:

Mathieu Morlighem

Message:

CHG: P0 controls now working in 3d

Location:

issm/trunk-jpl/src/c

Files:

• analyses/AdjointHorizAnalysis.cpp (modified) (12 diffs)
• classes/Elements/Penta.cpp (modified) (1 diff)
• classes/Inputs/PentaInput.cpp (modified) (1 diff)

issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp

@@ -1187 +1187 @@
         Element* basalelement;
 
+        if(control_interp==P0Enum) _error_("cannot require regularization for P0 controls");
+
         /*Get basal element*/
         element->FindParam(&domaintype,DomainTypeEnum);
@@ -1375 +1377 @@
         /*Intermediaries*/
         int      domaintype;
+
+        if(control_interp==P0Enum) _error_("cannot require regularization for P0 controls");
 
         /*Get basal element*/
@@ -1440 +1444 @@
         /*Intermediaries*/
         int      domaintype;
+        if(control_interp!=P1Enum) _error_("not implemented yet...");
 
         /*Get basal element*/
@@ -1502 +1507 @@
         /*Intermediaries*/
         int      domaintype,dim;
+        if(control_interp!=P1Enum) _error_("not implemented yet...");
 
         /*Get domaintype*/
@@ -1579 +1585 @@
         int      domaintype,dim;
         Element* basalelement;
+        if(control_interp!=P1Enum) _error_("not implemented yet...");
 
         /*Get basal element*/
@@ -1661 +1668 @@
         /*return if floating (gradient is 0)*/
         if(element->IsAllFloating()) return;
+        if(control_interp==P0Enum) _error_("cannot require regularization for P0 controls");
 
         /*Intermediaries*/
@@ -1755 +1763 @@
         if(element->IsAllFloating()) return;
         if(!element->IsOnBase()) return;
+        if(control_interp!=P1Enum) _error_("not implemented yet...");
 
         /*Intermediaries*/
@@ -1934 +1943 @@
 
                 /*Build gradient vector (actually -dJ/dD): */
-                for(int i=0;i<numvertices;i++){
-                        if(domaintype!=Domain2DverticalEnum) ge[i]+=-2.*drag*dalpha2dk*((lambda*vx+mu*vy))*Jdet*gauss->weight*basis[i];
-                        else ge[i]+=-2.*drag*dalpha2dk*(lambda*vx)*Jdet*gauss->weight*basis[i];
-                        _assert_(!xIsNan<IssmDouble>(ge[i]));
-                }
-        }
-        gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
+                if(control_interp==P1Enum){
+                        for(int i=0;i<numvertices;i++){
+                                if(domaintype!=Domain2DverticalEnum) ge[i]+=-2.*drag*dalpha2dk*((lambda*vx+mu*vy))*Jdet*gauss->weight*basis[i];
+                                else ge[i]+=-2.*drag*dalpha2dk*(lambda*vx)*Jdet*gauss->weight*basis[i];
+                                _assert_(!xIsNan<IssmDouble>(ge[i]));
+                        }
+                }
+                else if(control_interp==P0Enum){
+                        if(domaintype!=Domain2DverticalEnum) ge[0]+=-2.*drag*dalpha2dk*((lambda*vx+mu*vy))*Jdet*gauss->weight;
+                        else ge[0]+=-2.*drag*dalpha2dk*(lambda*vx)*Jdet*gauss->weight;
+                        _assert_(!xIsNan<IssmDouble>(ge[0]));
+                }
+                else{
+                        _error_("not supported");
+                }
+        }
+        if(control_interp==P1Enum){
+                gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
+        }
+        else if(control_interp==P0Enum){
+                gradient->SetValue(vertexpidlist[0],ge[0],ADD_VAL);
+        }
+        else{
+                _error_("not supported");
+        }
 
         /*Clean up and return*/
@@ -2073 +2100 @@
         if(element->IsAllFloating()) return;
         if(!element->IsOnBase()) return;
+        if(control_interp!=P1Enum) _error_("not implemented yet...");
 
         /*Intermediaries*/
@@ -2171 +2199 @@
         if(element->IsAllFloating()) return;
         if(!element->IsOnBase()) return;
+        if(control_interp!=P1Enum) _error_("not implemented yet...");
 
         /*Intermediaries*/
@@ -2241 +2270 @@
         /*return if floating (gradient is 0)*/
         if(element->IsAllFloating()) return;
+        if(control_interp!=P1Enum) _error_("not implemented yet...");
 
         /*Intermediaries*/
@@ -2366 +2396 @@
         int      domaintype,dim;
         Element* basalelement;
+        if(control_interp!=P1Enum) _error_("not implemented yet...");
 
         /*Get basal element*/

issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

@@ -220 +220 @@
         /*Call inputs method*/
         switch(interpolation_enum){
+                case P0Enum:
+                        control_input->SetControl(interpolation_enum,1,&this->lid,values,values_min,values_max);
+                        break;
                 case P1Enum:
                         control_input->SetControl(interpolation_enum,NUMVERTICES,&vertexlids[0],values,values_min,values_max);

issm/trunk-jpl/src/c/classes/Inputs/PentaInput.cpp

@@ -177 +177 @@
                 }
         }
+        else if(interp_in==P0Enum && this->interpolation==P0Enum){
+                _assert_(this->N==1);
+                for(int i=0;i<numindices;i++){
+                        int row = indices[i];
+                        _assert_(row>=0);
+                        _assert_(row<this->M);
+                        this->values[row] = values_in[i];
+                }
+        }
         else if(this->interpolation!=P1Enum && interp_in==P1Enum){
                 this->Reset(interp_in);