Changeset 13721 for issm/trunk-jpl/src/c/classes/FemModel.cpp

Timestamp:

10/17/12 15:59:20 (12 years ago)

Author:

Mathieu Morlighem

Message:

CHG: moved MassFluxx from modules to FemModel method

File:

• issm/trunk-jpl/src/c/classes/FemModel.cpp (modified) (6 diffs)

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

@@ -228 +228 @@
 }
 /*}}}*/
+/*FUNCTION FemModel::SetCurrentConfiguration(int configuration_type,int analysis_type){{{*/
+void FemModel::SetCurrentConfiguration(int configuration_type,int analysis_type){
+
+        /*Use configuration_type to setup the analysis counter, the configurations of objects etc ... but use 
+         * analysis_type to drive the element numerics. This allows for use of 1 configuration_type for several 
+         * analyses. For example: do a SurfaceSlopeX, SurfaceSlopeY, BedSlopeX and BedSlopeY analysis using the 
+         * Slope configuration.*/
+
+        int found=-1;
+        for(int i=0;i<nummodels;i++){
+                if (analysis_type_list[i]==configuration_type){
+                        found=i;
+                        break;
+                }
+        }
+        if(found!=-1) analysis_counter=found;
+        else _error_("Could not find alias for analysis_type " << EnumToStringx(configuration_type) << " in list of FemModel analyses");
+
+        /*Now, plug analysis_counter and analysis_type inside the parameters: */
+        this->parameters->SetParam(analysis_counter,AnalysisCounterEnum);
+        this->parameters->SetParam(analysis_type,AnalysisTypeEnum);
+        this->parameters->SetParam(configuration_type,ConfigurationTypeEnum);
+
+        /*configure elements, loads and nodes, for this new analysis: */
+        this->elements->SetCurrentConfiguration(elements,loads, nodes,vertices, materials,parameters);
+        this->nodes->SetCurrentConfiguration(elements,loads, nodes,vertices, materials,parameters);
+        this->loads->SetCurrentConfiguration(elements, loads, nodes,vertices, materials,parameters);
+
+        #ifdef _HAVE_PETSC_
+        /*take care of petsc options, that depend on this analysis type (present only after model processor)*/
+        if(this->parameters->Exist(PetscOptionsStringsEnum)){
+                PetscOptionsFromAnalysis(this->parameters,analysis_type);
+                if(VerboseSolver()) _pprintLine_("      petsc Options set for analysis type: " << EnumToStringx(analysis_type));
+        }
+        #endif
+
+}
+/*}}}*/
+/*FUNCTION FemModel::SetCurrentConfiguration(int configuration_type){{{*/
+void FemModel::SetCurrentConfiguration(int configuration_type){
+        this->SetCurrentConfiguration(configuration_type,configuration_type);
+}
+/*}}}*/
 /*FUNCTION FemModel::Solve {{{*/
 void FemModel::Solve(void){
@@ -280 +323 @@
 /*}}}*/
 
-/*Numerics: */
-/*FUNCTION FemModel::SetCurrentConfiguration(int configuration_type,int analysis_type){{{*/
-void FemModel::SetCurrentConfiguration(int configuration_type,int analysis_type){
-
-        /*Use configuration_type to setup the analysis counter, the configurations of objects etc ... but use 
-         * analysis_type to drive the element numerics. This allows for use of 1 configuration_type for several 
-         * analyses. For example: do a SurfaceSlopeX, SurfaceSlopeY, BedSlopeX and BedSlopeY analysis using the 
-         * Slope configuration.*/
-
-        int found=-1;
-        for(int i=0;i<nummodels;i++){
-                if (analysis_type_list[i]==configuration_type){
-                        found=i;
-                        break;
-                }
-        }
-        if(found!=-1) analysis_counter=found;
-        else _error_("Could not find alias for analysis_type " << EnumToStringx(configuration_type) << " in list of FemModel analyses");
-
-        /*Now, plug analysis_counter and analysis_type inside the parameters: */
-        this->parameters->SetParam(analysis_counter,AnalysisCounterEnum);
-        this->parameters->SetParam(analysis_type,AnalysisTypeEnum);
-        this->parameters->SetParam(configuration_type,ConfigurationTypeEnum);
-
-        /*configure elements, loads and nodes, for this new analysis: */
-        this->elements->SetCurrentConfiguration(elements,loads, nodes,vertices, materials,parameters);
-        this->nodes->SetCurrentConfiguration(elements,loads, nodes,vertices, materials,parameters);
-        this->loads->SetCurrentConfiguration(elements, loads, nodes,vertices, materials,parameters);
-
-        #ifdef _HAVE_PETSC_
-        /*take care of petsc options, that depend on this analysis type (present only after model processor)*/
-        if(this->parameters->Exist(PetscOptionsStringsEnum)){
-                PetscOptionsFromAnalysis(this->parameters,analysis_type);
-                if(VerboseSolver()) _pprintLine_("      petsc Options set for analysis type: " << EnumToStringx(analysis_type));
-        }
-        #endif
-
-}
-/*}}}*/
-/*FUNCTION FemModel::SetCurrentConfiguration(int configuration_type){{{*/
-void FemModel::SetCurrentConfiguration(int configuration_type){
-        this->SetCurrentConfiguration(configuration_type,configuration_type);
-}
-/*}}}*/
-
 /*Modules:*/
-int FemModel::UpdateVertexPositionsx(void){ /*{{{*/
+int  FemModel::UpdateVertexPositionsx(void){ /*{{{*/
 
         int     i;
@@ -378 +376 @@
         NodesDofx(nodes,parameters,analysis_type);
 
-}
-/*}}}*/
-void FemModel::TimeAdaptx(IssmDouble* pdt){/*{{{*/
-
-        int      i;
-
-        /*output: */
-        IssmDouble   dt;
-
-        /*intermediary: */
-        Element *element     = NULL;
-        IssmDouble   min_dt      = 0;
-        IssmDouble   node_min_dt = 0;
-
-        /*Go through elements, and figure out the minimum of the time steps for each element (using CFL criterion): */
-        element=(Element*)elements->GetObjectByOffset(0); min_dt=element->TimeAdapt();
-
-        for (i=1;i<elements->Size();i++){
-                element=(Element*)elements->GetObjectByOffset(i);
-                dt=element->TimeAdapt();
-                if(dt<min_dt)min_dt=dt;
-        }
-
-        /*Figure out minimum across the cluster: */
-        #ifdef _HAVE_MPI_
-        MPI_Reduce (&min_dt,&node_min_dt,1,MPI_DOUBLE,MPI_MIN,0,IssmComm::GetComm() );
-        MPI_Bcast(&node_min_dt,1,MPI_DOUBLE,0,IssmComm::GetComm());
-        min_dt=node_min_dt;
-        #endif
-
-        /*Assign output pointers:*/
-        *pdt=min_dt;
 }
 /*}}}*/
@@ -440 +406 @@
                 case MaxVzEnum:                  MaxVzx(                   responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
                 case MaxAbsVzEnum:               MaxAbsVzx(                responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
-                case MassFluxEnum:               MassFluxx(                responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;
+                case MassFluxEnum:               this->MassFluxx(          responses,process_units); break;
                 case SurfaceAbsVelMisfitEnum:    SurfaceAbsVelMisfitx(     responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break;
                 case SurfaceRelVelMisfitEnum:    SurfaceRelVelMisfitx(     responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break;
@@ -543 +509 @@
 }
 /*}}}*/
+void FemModel::TimeAdaptx(IssmDouble* pdt){/*{{{*/
+
+        int      i;
+
+        /*output: */
+        IssmDouble   dt;
+
+        /*intermediary: */
+        Element *element     = NULL;
+        IssmDouble   min_dt      = 0;
+        IssmDouble   node_min_dt = 0;
+
+        /*Go through elements, and figure out the minimum of the time steps for each element (using CFL criterion): */
+        element=(Element*)elements->GetObjectByOffset(0); min_dt=element->TimeAdapt();
+
+        for (i=1;i<elements->Size();i++){
+                element=(Element*)elements->GetObjectByOffset(i);
+                dt=element->TimeAdapt();
+                if(dt<min_dt)min_dt=dt;
+        }
+
+        /*Figure out minimum across the cluster: */
+#ifdef _HAVE_MPI_
+        MPI_Reduce (&min_dt,&node_min_dt,1,MPI_DOUBLE,MPI_MIN,0,IssmComm::GetComm() );
+        MPI_Bcast(&node_min_dt,1,MPI_DOUBLE,0,IssmComm::GetComm());
+        min_dt=node_min_dt;
+#endif
+
+        /*Assign output pointers:*/
+        *pdt=min_dt;
+}
+/*}}}*/
+#ifdef _HAVE_RESPONSES_
+void FemModel::MassFluxx(IssmDouble* pmass_flux,bool process_units){/*{{{*/
+
+        int          i,j;
+        Element     *element       = NULL;
+        int          element_id;
+        bool         ispresent     = false;
+        IssmDouble   mass_flux     = 0;
+        IssmDouble   all_mass_flux = 0;
+        int          counter;
+        IssmDouble **array         = NULL;
+        int          M;
+        int         *mdims_array   = NULL;
+        int         *ndims_array   = NULL;
+        IssmDouble  *segments      = NULL;
+        int          num_segments;
+
+        /*First, figure out which segment to compute our mass flux on. Start with retrieving qmu_mass_flux_segments: */
+        this->parameters->FindParam(&ispresent,MassFluxSegmentsPresentEnum);
+        if(!ispresent)_error_("no mass flux segments available!");
+        this->parameters->FindParam(&array,&M,&mdims_array,&ndims_array,MassFluxSegmentsEnum);
+
+        /*Retrieve index of segments being used for MassFlux computation: */
+        parameters->FindParam(&counter,IndexEnum);
+
+        /*retrieve segments from array: */
+        segments     = array[counter-1]; //matlab to "C" indexing
+        num_segments = mdims_array[counter-1];
+
+        /*Go through segments, and then elements, and figure out which elements belong to a segment. 
+         * When we find one, use the element to compute the mass flux on the segment: */
+        for(i=0;i<num_segments;i++){
+                element_id=reCast<int,IssmDouble>(*(segments+5*i+4));
+                for(j=0;j<elements->Size();j++){
+                        element=(Element*)this->elements->GetObjectByOffset(j);
+                        if (element->Id()==element_id){
+                                /*We found the element which owns this segment, use it to compute the mass flux: */
+                                mass_flux+=element->MassFlux(segments+5*i+0,process_units);
+                                break;
+                        }
+                }
+        }
+
+#ifdef _HAVE_MPI_
+        MPI_Allreduce ( (void*)&mass_flux,(void*)&all_mass_flux,1,MPI_DOUBLE,MPI_SUM,IssmComm::GetComm());
+        mass_flux=all_mass_flux;
+#endif
+
+        /*Free ressources:*/
+        for(i=0;j<M;i++){
+                IssmDouble* matrix=array[j];
+                xDelete<IssmDouble>(matrix);
+        }
+        xDelete<int>(mdims_array);
+        xDelete<int>(ndims_array);
+        xDelete<IssmDouble*>(array);
+
+        /*Assign output pointers: */
+        *pmass_flux=mass_flux;
+
+}/*}}}*/
+#endif
 #ifdef _HAVE_CONTROL_
 void FemModel::ThicknessAbsGradientx( IssmDouble* pJ, bool process_units, int weight_index){/*{{{*/
@@ -601 +661 @@
 /*}}}*/
 #endif
-
 #ifdef  _HAVE_DAKOTA_
 void FemModel::DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses){/*{{{*/