Changeset 24490

Timestamp:

12/20/19 14:32:04 (5 years ago)

Author:

tsantos

Message:

CHG: added thickness average method for Discontinuous Galerkin (stabilization in Mass Transport Analysis). Not totally tested, just keeping the code

Location:

issm/trunk-jpl/src/c

Files:

• classes/FemModel.cpp (modified) (1 diff)
• classes/FemModel.h (modified) (1 diff)
• cores/masstransport_core.cpp (modified) (1 diff)

issm/trunk-jpl/src/c/classes/FemModel.cpp

@@ -2761 +2761 @@
 }
 /*}}}*/
+void FemModel::ThicknessAverage(){/*{{{*/
+   
+        int elementswidth                   = this->GetElementsWidth();//just 2D mesh, tria elements
+   int numberofvertices                = this->vertices->NumberOfVertices();//total number of vertices
+
+   IssmDouble weight                   = 0.;
+   IssmDouble* totalweight             = NULL;
+        IssmDouble* Hserial                                             = NULL;
+   IssmDouble* H                       = xNew<IssmDouble>(elementswidth);
+   Vector<IssmDouble>* vecH                             = new Vector<IssmDouble>(numberofvertices);
+   Vector<IssmDouble>* vectotalweight  = new Vector<IssmDouble>(numberofvertices);
+
+   for(int i=0;i<this->elements->Size();i++){
+      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
+     
+                /*check if there is ice in this element*/
+                if(!element->IsIceInElement()) continue;  
+        
+                /*get H on the vertices*/
+                element->GetInputListOnVertices(H,ThicknessEnum);
+
+      /*weight to calculate the smoothed H*/
+      weight=1.;//simple average
+
+                /*add in the serial vector*/
+      vecH->SetValue(element->vertices[0]->Sid(),weight*H[0],ADD_VAL);
+      vecH->SetValue(element->vertices[1]->Sid(),weight*H[1],ADD_VAL);
+      vecH->SetValue(element->vertices[2]->Sid(),weight*H[2],ADD_VAL);
+      /*total weight*/
+      vectotalweight->SetValue(element->vertices[0]->Sid(),weight,ADD_VAL);
+      vectotalweight->SetValue(element->vertices[1]->Sid(),weight,ADD_VAL);
+      vectotalweight->SetValue(element->vertices[2]->Sid(),weight,ADD_VAL);
+   }
+
+   /*Assemble and serialize*/
+   vecH->Assemble();
+   vectotalweight->Assemble();
+   Hserial=vecH->ToMPISerial();
+   totalweight=vectotalweight->ToMPISerial();
+
+   /*Divide for the total weight*/
+   for(int i=0;i<numberofvertices;i++){
+      _assert_(totalweight[i]>0);
+      Hserial[i]=Hserial[i]/totalweight[i];
+   }
+
+   /*Set element inputs*/
+   for(int i=0;i<this->elements->Size();i++){
+      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
+                H[0]=Hserial[element->vertices[0]->Sid()];
+                H[1]=Hserial[element->vertices[1]->Sid()];
+                H[2]=Hserial[element->vertices[2]->Sid()];
+                element->AddInput2(ThicknessEnum,H,P1Enum);
+        }
+
+        /*Cleanup*/
+   delete vecH;
+   delete vectotalweight;
+   xDelete<IssmDouble>(H);
+   xDelete<IssmDouble>(Hserial);
+   xDelete<IssmDouble>(totalweight);
+}
+/*}}}*/
 void FemModel::ThicknessPositivex(IssmDouble* pJ){/*{{{*/
 

issm/trunk-jpl/src/c/classes/FemModel.h

@@ -151 +151 @@
                 void OmegaAbsGradientx( IssmDouble* pJ);
                 void EtaDiffx( IssmDouble* pJ);
+                void ThicknessAverage();
                 void ThicknessAbsGradientx( IssmDouble* pJ);
                 void ThicknessPositivex(IssmDouble* pJ);

issm/trunk-jpl/src/c/cores/masstransport_core.cpp

@@ -70 +70 @@
                 else{
                         solutionsequence_linear(femmodel);
+                        // ThicknessAverage: method not totally tested
+                        //if(stabilization==3){
+                        //      if(VerboseSolution()) _printf0_("   call thickness average core\n");
+                        //      femmodel->ThicknessAverage();
+                        //}
                 }
                 femmodel->parameters->SetParam(ThicknessEnum,InputToExtrudeEnum);