Changeset 18058 for issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp

Timestamp:

05/25/14 18:01:37 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: working on moving gradient to adjoint

File:

• issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (modified) (1 diff)

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

@@ -974 +974 @@
 }/*}}}*/
 void AdjointHorizAnalysis::GradientJDragGradient(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
-        _error_("not implemented yet");
+
+        /*return if floating*/
+        if(element->IsFloating())return;
+
+        /*Intermediaries*/
+        int      domaintype,dim;
+        Element* basalelement;
+
+        /*Get basal element*/
+        element->FindParam(&domaintype,DomainTypeEnum);
+        switch(domaintype){
+                case Domain2DhorizontalEnum:
+                        basalelement = element;
+                        dim          = 2;
+                        break;
+                case Domain2DverticalEnum:
+                        if(!element->IsOnBase()) return;
+                        basalelement = element->SpawnBasalElement();
+                        dim          = 1;
+                        break;
+                case Domain3DEnum:
+                        if(!element->IsOnBase()) return;
+                        basalelement = element->SpawnBasalElement();
+                        dim          = 2;
+                        break;
+                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
+        }
+
+        /*Intermediaries*/
+        IssmDouble Jdet,weight;
+        IssmDouble dk[3]; 
+        IssmDouble *xyz_list= NULL;
+
+        /*Fetch number of vertices for this finite element*/
+        int numvertices = basalelement->GetNumberOfVertices();
+
+        /*Initialize some vectors*/
+        IssmDouble* dbasis        = xNew<IssmDouble>(2*numvertices);
+        IssmDouble* ge            = xNew<IssmDouble>(numvertices);
+        int*        vertexpidlist = xNew<int>(numvertices);
+
+        /*Retrieve all inputs we will be needing: */
+        basalelement->GetVerticesCoordinates(&xyz_list);
+        basalelement->GradientIndexing(&vertexpidlist[0],control_index);
+        Input* dragcoefficient_input = basalelement->GetInput(FrictionCoefficientEnum);                _assert_(dragcoefficient_input);
+        Input* weights_input         = basalelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=basalelement->NewGauss(2);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+
+                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                basalelement->NodalFunctionsP1Derivatives(dbasis,xyz_list,gauss);
+                weights_input->GetInputValue(&weight,gauss,DragCoefficientAbsGradientEnum);
+
+                /*Build alpha_complement_list: */
+                dragcoefficient_input->GetInputDerivativeValue(&dk[0],xyz_list,gauss);
+
+                /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
+                for(int i=0;i<numvertices;i++){
+                        if(dim==2){
+                                ge[i]+=-weight*Jdet*gauss->weight*(dbasis[0*numvertices+i]*dk[0]+dbasis[1*numvertices+i]*dk[1]);
+                        }
+                        else{
+                                ge[i]+=-weight*Jdet*gauss->weight*dbasis[0*numvertices+i]*dk[0];
+                        }
+                        _assert_(!xIsNan<IssmDouble>(ge[i]));
+                }
+        }
+        gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(xyz_list);
+        xDelete<IssmDouble>(dbasis);
+        xDelete<IssmDouble>(ge);
+        xDelete<int>(vertexpidlist);
+        delete gauss;
+        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
+
 }/*}}}*/
 void AdjointHorizAnalysis::GradientJBGradient(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/