[Core][MeshingApplication] Extend gradient process to other geometries. Reduce the number of template arguments to reduce compilation time

applications/MeshingApplication/custom_processes/metrics_hessian_process.h

@@ -53,7 +53,6 @@ namespace Kratos
  * @brief This class is can be used to compute the metrics of the model part with an Hessian approach
  * @author Vicente Mataix Ferrandiz
  */
-template<SizeType TDim, class TVarType>
 class KRATOS_API(MESHING_APPLICATION) ComputeHessianSolMetricProcess
     : public Process
 {
@@ -73,12 +72,6 @@ class KRATOS_API(MESHING_APPLICATION) ComputeHessianSolMetricProcess
     /// The index type definition
     typedef std::size_t                                                            IndexType;
 
-    /// The type of array considered for the tensor
-    typedef typename std::conditional<TDim == 2, array_1d<double, 3>, array_1d<double, 6>>::type TensorArrayType;
-
-    /// Matrix type definition
-    typedef BoundedMatrix<double, TDim, TDim> MatrixType;
-
     /// Pointer definition of ComputeHessianSolMetricProcess
     KRATOS_CLASS_POINTER_DEFINITION(ComputeHessianSolMetricProcess);
 
@@ -98,15 +91,26 @@ class KRATOS_API(MESHING_APPLICATION) ComputeHessianSolMetricProcess
     // Constructor
 
     /**
-     * @brief This is the default constructor
+     * @brief This is the default constructor (double)
      * @param rThisModelPart The model part to be computed
      * @param rVariable The variable to compute
      * @param ThisParameters The input parameters
      */
+    ComputeHessianSolMetricProcess(
+        ModelPart& rThisModelPart,
+        Variable<double>& rVariable,
+        Parameters ThisParameters = Parameters(R"({})")
+        );
 
+    /**
+     * @brief This is the default constructor (component)
+     * @param rThisModelPart The model part to be computed
+     * @param rVariable The variable to compute
+     * @param ThisParameters The input parameters
+     */
     ComputeHessianSolMetricProcess(
         ModelPart& rThisModelPart,
-        TVarType& rVariable,
+        ComponentType& rVariable,
         Parameters ThisParameters = Parameters(R"({})")
         );
 
@@ -208,18 +212,21 @@ class KRATOS_API(MESHING_APPLICATION) ComputeHessianSolMetricProcess
     ///@name Private member Variables
     ///@{
 
-    ModelPart& mThisModelPart;                           /// The model part to compute
-    TVarType& mVariable;                                 /// The variable to calculate the hessian
-    std::string mRatioReferenceVariable = "DISTANCE";    /// Variable used to compute the anisotropic ratio
-    double mMinSize;                                     /// The minimal size of the elements
-    double mMaxSize;                                     /// The maximal size of the elements
-    bool mEnforceCurrent;                                /// With this we choose if we inforce the current nodal size (NODAL_H)
-    bool mEstimateInterpError;                           /// If the error of interpolation will be estimated
-    double mInterpError;                                 /// The error of interpolation allowed
-    double mMeshConstant;                                /// The mesh constant to remesh (depends of the element type)
-    double mAnisotropicRatio;                            /// The minimal anisotropic ratio (0 < ratio < 1)
-    double mBoundLayer;                                  /// The boundary layer limit distance
-    Interpolation mInterpolation;                        /// The interpolation type
+    ModelPart& mThisModelPart;                                  /// The model part to compute
+    std::vector<Variable<double>*> mrOriginVariableDoubleList;  /// The scalar variable list to compute
+    std::vector<ComponentType*> mrOriginVariableComponentsList; /// The scalar variable list to compute (components)
+
+    // TODO: Replace for Parameters
+    std::string mRatioReferenceVariable = "DISTANCE";          /// Variable used to compute the anisotropic ratio
+    double mMinSize;                                           /// The minimal size of the elements
+    double mMaxSize;                                           /// The maximal size of the elements
+    bool mEnforceCurrent;                                      /// With this we choose if we inforce the current nodal size (NODAL_H)
+    bool mEstimateInterpError;                                 /// If the error of interpolation will be estimated
+    double mInterpError;                                       /// The error of interpolation allowed
+    double mMeshConstant;                                      /// The mesh constant to remesh (depends of the element type)
+    double mAnisotropicRatio;                                  /// The minimal anisotropic ratio (0 < ratio < 1)
+    double mBoundLayer;                                        /// The boundary layer limit distance
+    Interpolation mInterpolation;                              /// The interpolation type
 
     ///@}
     ///@name Private Operators
@@ -236,12 +243,15 @@ class KRATOS_API(MESHING_APPLICATION) ComputeHessianSolMetricProcess
      * @param AnisotropicRatio The anisotropic ratio
      * @param ElementMinSize The min size of element
      * @param ElementMaxSize The maximal size of the elements
+     * @param NodalH The size of the local node
      */
+    template<SizeType TDim>
     array_1d<double, 3 * (TDim - 1)> ComputeHessianMetricTensor(
         const Vector& rHessian,
         const double AnisotropicRatio,
         const double ElementMinSize, // This way we can impose as minimum as the previous size if we desire
-        const double ElementMaxSize // This way we can impose as maximum as the previous size if we desire
+        const double ElementMaxSize, // This way we can impose as maximum as the previous size if we desire
+        const double NodalH
         );
 
     /**
@@ -280,6 +290,22 @@ class KRATOS_API(MESHING_APPLICATION) ComputeHessianSolMetricProcess
         const Interpolation rInterpolation
         );
 
+    /**
+     * @brief This method is the responsible to compute the metric of the problem
+     */
+    template<SizeType TDim>
+    void CalculateMetric();
+
+    /**
+     * @brief This method provides the defaults parameters to avoid conflicts between the different constructors
+     */
+    Parameters GetDefaultParameters() const;
+
+    /**
+     * @brief This method provides the defaults parameters to avoid conflicts between the different constructors
+     */
+    void InitializeVariables(Parameters ThisParameters);
+
     ///@}
     ///@name Private  Access
     ///@{
@@ -316,14 +342,12 @@ class KRATOS_API(MESHING_APPLICATION) ComputeHessianSolMetricProcess
 ///@{
 
 /// input stream function
-template<unsigned int TDim, class TVarType>
 inline std::istream& operator >> (std::istream& rIStream,
-                                  ComputeHessianSolMetricProcess<TDim, TVarType>& rThis);
+                                  ComputeHessianSolMetricProcess& rThis);
 
 /// output stream function
-template<unsigned int TDim, class TVarType>
 inline std::ostream& operator << (std::ostream& rOStream,
-                                  const ComputeHessianSolMetricProcess<TDim, TVarType>& rThis)
+                                  const ComputeHessianSolMetricProcess& rThis)
 {
     rThis.PrintInfo(rOStream);
     rOStream << std::endl;

applications/MeshingApplication/custom_python/add_processes_to_python.cpp

@@ -89,27 +89,18 @@ void  AddProcessesToPython(pybind11::module& m)
     .def(py::init<ModelPart&, const Variable<array_1d<double,3>>, Parameters>())
     ;
 
-    // HESSIAN DOUBLE
-    py::class_<ComputeHessianSolMetricProcess<2, Variable<double>>, ComputeHessianSolMetricProcess<2, Variable<double>>::Pointer, Process>(m, "ComputeHessianSolMetricProcess2D")
+    // HESSIAN PROCESS
+    py::class_<ComputeHessianSolMetricProcess, ComputeHessianSolMetricProcess::Pointer, Process>(m, "ComputeHessianSolMetricProcess")
     .def(py::init<ModelPart&, Variable<double>&>())
     .def(py::init<ModelPart&, Variable<double>&, Parameters>())
-    ;
-
-    py::class_<ComputeHessianSolMetricProcess<3, Variable<double>>, ComputeHessianSolMetricProcess<3, Variable<double>>::Pointer, Process>(m, "ComputeHessianSolMetricProcess3D")
-    .def(py::init<ModelPart&, Variable<double>&>())
-    .def(py::init<ModelPart&, Variable<double>&, Parameters>())
-    ;
-
-    // HESSIAN ARRAY 1D
-    py::class_<ComputeHessianSolMetricProcess<2, ComponentType>, ComputeHessianSolMetricProcess<2, ComponentType>::Pointer, Process>(m, "ComputeHessianSolMetricProcessComp2D")
     .def(py::init<ModelPart&, ComponentType&>())
     .def(py::init<ModelPart&, ComponentType&, Parameters>())
     ;
 
-    py::class_<ComputeHessianSolMetricProcess<3, ComponentType>, ComputeHessianSolMetricProcess<3, ComponentType>::Pointer, Process>(m, "ComputeHessianSolMetricProcessComp3D")
-    .def(py::init<ModelPart&, ComponentType&>())
-    .def(py::init<ModelPart&, ComponentType&, Parameters>())
-    ;
+    m.attr("ComputeHessianSolMetricProcess2D") = m.attr("ComputeHessianSolMetricProcess");
+    m.attr("ComputeHessianSolMetricProcess3D") = m.attr("ComputeHessianSolMetricProcess");
+    m.attr("ComputeHessianSolMetricProcessComp2D") = m.attr("ComputeHessianSolMetricProcess");
+    m.attr("ComputeHessianSolMetricProcessComp3D") = m.attr("ComputeHessianSolMetricProcess");
 
     // ERROR
     py::class_<MetricErrorProcess<2>, MetricErrorProcess<2>::Pointer, Process>(m, "MetricErrorProcess2D")

applications/MeshingApplication/python_scripts/mmg_process.py

@@ -123,12 +123,21 @@ def __init__(self, Model, settings ):
         }
         """)
 
+        # Identify the dimension first
+        if not settings.Has("model_part_name"):
+            settings.AddValue("model_part_name", default_parameters["model_part_name"])
+        self.dim = Model[settings["model_part_name"].GetString()].ProcessInfo[KratosMultiphysics.DOMAIN_SIZE]
+        # The mesh dependent constant depends on dimension
+        if self.dim == 2:
+            default_parameters["hessian_strategy_parameters"]["mesh_dependent_constant"].SetDouble(2.0/9.0)
+        else:
+            default_parameters["hessian_strategy_parameters"]["mesh_dependent_constant"].SetDouble(9.0/32.0)
+
         # Overwrite the default settings with user-provided parameters
         self.settings = settings
         self.settings.RecursivelyValidateAndAssignDefaults(default_parameters)
 
         self.model_part= Model[self.settings["model_part_name"].GetString()]
-        self.dim = self.model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE]
 
         self.enforce_current = self.settings["enforce_current"].GetBool()
 

applications/MeshingApplication/tests/cpp_tests/test_metric_process.cpp

@@ -244,8 +244,9 @@ namespace Kratos
             this_model_part.AddNodalSolutionStepVariable(DISTANCE_GRADIENT);
 
             auto& process_info = this_model_part.GetProcessInfo();
-            process_info[STEP] = 1;
-            process_info[NL_ITERATION_NUMBER] = 1;
+            process_info.SetValue(DOMAIN_SIZE, 2);
+            process_info.SetValue(STEP, 1);
+            process_info.SetValue(NL_ITERATION_NUMBER, 1);
 
             Create2DGeometry(this_model_part, "Element2D3N");
 
@@ -254,10 +255,11 @@ namespace Kratos
                 auto it_node = this_model_part.Nodes().begin() + i_node;
                 it_node->FastGetSolutionStepValue(DISTANCE) = (it_node->X() == 1.0) ? 0.0 : 1.0;
                 it_node->SetValue(NODAL_H, 1.0);
+                it_node->SetValue(NODAL_AREA, 0.0);
                 it_node->SetValue(METRIC_TENSOR_2D, ZeroVector(3));
             }
 
-            typedef ComputeNodalGradientProcess<2, Variable<double>, Historical> GradientType;
+            typedef ComputeNodalGradientProcess<ComputeNodalGradientProcessSettings::SaveAsHistoricalVariable> GradientType;
             GradientType gradient_process = GradientType(this_model_part, DISTANCE, DISTANCE_GRADIENT, NODAL_AREA);
             gradient_process.Execute();
 
@@ -293,8 +295,9 @@ namespace Kratos
             this_model_part.AddNodalSolutionStepVariable(DISTANCE_GRADIENT);
 
             auto& process_info = this_model_part.GetProcessInfo();
-            process_info[STEP] = 1;
-            process_info[NL_ITERATION_NUMBER] = 1;
+            process_info.SetValue(DOMAIN_SIZE, 3);
+            process_info.SetValue(STEP, 1);
+            process_info.SetValue(NL_ITERATION_NUMBER, 1);
 
             Create3DGeometry(this_model_part, "Element3D4N");
 
@@ -303,11 +306,12 @@ namespace Kratos
                 auto it_node = this_model_part.Nodes().begin() + i_node;
                 it_node->FastGetSolutionStepValue(DISTANCE) = (it_node->X() == 1.0) ? 0.0 : 1.0;
                 it_node->SetValue(NODAL_H, 1.0);
+                it_node->SetValue(NODAL_AREA, 0.0);
                 it_node->SetValue(METRIC_TENSOR_3D, ZeroVector(6));
             }
 
             // Compute gradient
-            typedef ComputeNodalGradientProcess<3, Variable<double>, Historical> GradientType;
+            typedef ComputeNodalGradientProcess<ComputeNodalGradientProcessSettings::SaveAsHistoricalVariable> GradientType;
             GradientType gradient_process = GradientType(this_model_part, DISTANCE, DISTANCE_GRADIENT, NODAL_AREA);
             gradient_process.Execute();
 
@@ -347,8 +351,9 @@ namespace Kratos
             this_model_part.AddNodalSolutionStepVariable(DISTANCE_GRADIENT);
 
             auto& process_info = this_model_part.GetProcessInfo();
-            process_info[STEP] = 1;
-            process_info[NL_ITERATION_NUMBER] = 1;
+            process_info.SetValue(DOMAIN_SIZE, 2);
+            process_info.SetValue(STEP, 1);
+            process_info.SetValue(NL_ITERATION_NUMBER, 1);
 
             Create2DGeometry(this_model_part, "Element2D3N");
 
@@ -361,7 +366,7 @@ namespace Kratos
             }
 
             // Compute metric
-            ComputeHessianSolMetricProcess<2, Variable<double>> hessian_process = ComputeHessianSolMetricProcess<2, Variable<double>>(this_model_part, DISTANCE);
+            auto hessian_process = ComputeHessianSolMetricProcess(this_model_part, DISTANCE);
             hessian_process.Execute();
 
 //             // DEBUG
@@ -392,8 +397,9 @@ namespace Kratos
             this_model_part.AddNodalSolutionStepVariable(DISTANCE_GRADIENT);
 
             auto& process_info = this_model_part.GetProcessInfo();
-            process_info[STEP] = 1;
-            process_info[NL_ITERATION_NUMBER] = 1;
+            process_info.SetValue(DOMAIN_SIZE, 3);
+            process_info.SetValue(STEP, 1);
+            process_info.SetValue(NL_ITERATION_NUMBER, 1);
 
             Create3DGeometry(this_model_part, "Element3D4N");
 
@@ -406,7 +412,7 @@ namespace Kratos
             }
 
             // Compute metric
-            ComputeHessianSolMetricProcess<3, Variable<double>> hessian_process = ComputeHessianSolMetricProcess<3, Variable<double>>(this_model_part, DISTANCE);
+            auto hessian_process = ComputeHessianSolMetricProcess(this_model_part, DISTANCE);
             hessian_process.Execute();
 
 //             // DEBUG
@@ -440,8 +446,9 @@ namespace Kratos
             this_model_part.AddNodalSolutionStepVariable(DISPLACEMENT);
 
             auto& process_info = this_model_part.GetProcessInfo();
-            process_info[STEP] = 1;
-            process_info[NL_ITERATION_NUMBER] = 1;
+            process_info.SetValue(DOMAIN_SIZE, 2);
+            process_info.SetValue(STEP, 1);
+            process_info.SetValue(NL_ITERATION_NUMBER, 1);
 
             // In case the StructuralMechanicsApplciation is not compiled we skip the test
             Properties::Pointer p_elem_prop = this_model_part.pGetProperties(0);
@@ -504,8 +511,9 @@ namespace Kratos
             this_model_part.AddNodalSolutionStepVariable(DISPLACEMENT);
 
             auto& process_info = this_model_part.GetProcessInfo();
-            process_info[STEP] = 1;
-            process_info[NL_ITERATION_NUMBER] = 1;
+            process_info.SetValue(DOMAIN_SIZE, 3);
+            process_info.SetValue(STEP, 1);
+            process_info.SetValue(NL_ITERATION_NUMBER, 1);
 
             // In case the StructuralMechanicsApplciation is not compiled we skip the test
             Properties::Pointer p_elem_prop = this_model_part.pGetProperties(0);

applications/MeshingApplication/tests/mmg_lagrangian_test/beam2D_line_load_test_parameters.json

@@ -87,16 +87,6 @@
         },
         "point_data_configuration"  : []
     },
-    "restart_options"          : {
-        "SaveRestart"      : false,
-        "RestartFrequency" : 0,
-        "LoadRestart"      : false,
-        "Restart_Step"     : 0
-    },
-    "constraints_data"         : {
-        "incremental_load"         : false,
-        "incremental_displacement" : false
-    },
     "recursive_remeshing_process" :[
     {
         "python_module"   : "mmg_process",

applications/MeshingApplication/tests/mmg_lagrangian_test/beam2D_line_load_test_result.sol

@@ -21,38 +21,38 @@ SolAtVertices
  1596.1682065273 0 1596.1682065273
  1596.1682065273 0 1596.1682065273
  1596.1682065273 0 1596.1682065273
- 1973.83483408483 0 1973.83483408483
+ 1741.2857324904 0 1741.2857324904
  1596.1682065273 0 1596.1682065273
  1596.1682065273 0 1596.1682065273
  1596.1682065273 0 1596.1682065273
  1596.1682065273 0 1596.1682065273
- 1783.02450156543 0 1783.02450156543
- 2233.44222914806 0 2233.44222914806
  1596.1682065273 0 1596.1682065273
+ 1888.41664848122 0 1888.41664848122
  1596.1682065273 0 1596.1682065273
- 2830.0529672361 0 2830.0529672361
- 2333.06689841094 0 2333.06689841094
- 2754.22625279034 0 2754.22625279034
- 3085.06572970157 0 3085.06572970157
- 2033.03549117724 0 2033.03549117724
- 3648.33349336322 0 3648.33349336322
- 3387.53606158952 0 3387.53606158952
- 3054.17230013612 0 3054.17230013612
- 3919.07288195414 0 3919.07288195414
- 4359.98062035213 0 4359.98062035213
- 4580.1152299552 0 4580.1152299552
- 4193.9182665218 0 4193.9182665218
- 4564.40173273184 0 4564.40173273184
- 5444.9092312109 0 5444.9092312109
- 6157.29896814932 0 6157.29896814932
- 6384.67282610918 0 6384.67282610918
- 5509.26950609488 0 5509.26950609488
- 5817.25316009115 0 5817.25316009115
- 6384.67282610918 0 6384.67282610918
- 6384.67282610918 0 6384.67282610918
- 6384.67282610918 0 6384.67282610918
- 6384.67282610918 0 6384.67282610918
+ 1596.1682065273 0 1596.1682065273
+ 2537.51413246724 0 2537.51413246724
+ 2010.22900468851 0 2010.22900468851
+ 2238.28410586729 0 2238.28410586729
+ 2763.41483935871 0 2763.41483935871
+ 1596.1682065273 0 1596.1682065273
+ 3299.94421406067 0 3299.94421406067
+ 2782.8761041239 0 2782.8761041239
+ 2324.43332491443 0 2324.43332491443
+ 3201.60086725174 0 3201.60086725174
+ 3500.57387518879 0 3500.57387518879
+ 4093.11138607511 0 4093.11138607511
+ 3173.77625263271 0 3173.77625263271
+ 3824.39707557243 0 3824.39707557243
+ 4223.97301372814 0 4223.97301372814
+ 4754.94642218925 0 4754.94642218925
+ 5513.76587138039 0 5513.76587138039
+ 4156.60334075396 0 4156.60334075396
+ 4969.4970415752 0 4969.4970415752
+ 5325.65916538894 0 5325.65916538894
+ 5933.73545574369 0 5933.73545574369
+ 5372.02361021606 0 5372.02361021606
  6384.67282610918 0 6384.67282610918
+ 6205.42714745037 0 6205.42714745037
  6384.67282610918 0 6384.67282610918
  6384.67282610918 0 6384.67282610918
  6384.67282610918 0 6384.67282610918