Merge branch 'master' into core/improve-document-test

README.md

@@ -1,6 +1,6 @@
 <p align=center><img height="72.125%" width="72.125%" src="https://raw.githubusercontent.com/KratosMultiphysics/Documentation/master/Wiki_files/Home/kratos.png"></p>
 
-[![License][license-image]][license] [![C++][c++-image]][c++standard] [![Github CI][Nightly-Build]][Nightly-link] [![DOI][DOI-image]][DOI] [![GitHub stars][stars-image]][stars] [![Twitter][twitter-image]][twitter]
+[![License][license-image]][license] [![C++][c++-image]][c++standard] [![Github CI][Nightly-Build]][Nightly-link] [![DOI][DOI-image]][DOI] [![GitHub stars][stars-image]][stars] [![Twitter][twitter-image]][twitter] [![Youtube][youtube-image]][youtube]
 
 [![Release][release-image]][releases]
 <a href="https://github.com/KratosMultiphysics/Kratos/releases/latest"><img src="https://img.shields.io/github/release-date/KratosMultiphysics/Kratos?label="></a>
@@ -10,7 +10,7 @@
 [![PyPI pyversions](https://img.shields.io/pypi/pyversions/KratosMultiphysics.svg)](https://pypi.org/project/KratosMultiphysics/)
 [![Downloads](https://pepy.tech/badge/KratosMultiphysics/month)](https://pepy.tech/project/KratosMultiphysics)
 
-[release-image]: https://img.shields.io/badge/release-9.3-green.svg?style=flat
+[release-image]: https://img.shields.io/badge/release-9.4-green.svg?style=flat
 [releases]: https://github.com/KratosMultiphysics/Kratos/releases
 
 [license-image]: https://img.shields.io/badge/license-BSD-green.svg?style=flat
@@ -31,6 +31,9 @@
 [twitter-image]: https://img.shields.io/twitter/follow/kratosmultiphys.svg?label=Follow&style=social
 [twitter]: https://twitter.com/kratosmultiphys
 
+[youtube-image]: https://badges.aleen42.com/src/youtube.svg
+[youtube]:https://www.youtube.com/@kratosmultiphysics3578
+
 _KRATOS Multiphysics_ ("Kratos") is a framework for building parallel, multi-disciplinary simulation software, aiming at modularity, extensibility, and high performance. Kratos is written in C++, and counts with an extensive Python interface. More in [Overview](https://github.com/KratosMultiphysics/Kratos/wiki/Overview)
 
 **Kratos** is **free** under BSD-4 [license](https://github.com/KratosMultiphysics/Kratos/wiki/Licence) and can be used even in comercial softwares as it is. Many of its main applications are also free and BSD-4 licensed but each derived application can have its own propietary license.

applications/ConvectionDiffusionApplication/convection_diffusion_application_variables.cpp

@@ -23,7 +23,6 @@ KRATOS_CREATE_VARIABLE(double, MEAN_SIZE)
 KRATOS_CREATE_VARIABLE(double, MEAN_VEL_OVER_ELEM_SIZE)
 KRATOS_CREATE_VARIABLE(double, MELT_TEMPERATURE_1)
 KRATOS_CREATE_VARIABLE(double, MELT_TEMPERATURE_2)
-KRATOS_CREATE_VARIABLE(double, PENALTY_DIRICHLET)
 KRATOS_CREATE_VARIABLE(double, EMBEDDED_SCALAR)
 KRATOS_CREATE_VARIABLE(double, PROJECTED_SCALAR1)
 KRATOS_CREATE_VARIABLE(double, TRANSFER_COEFFICIENT)

applications/ConvectionDiffusionApplication/convection_diffusion_application_variables.h

@@ -34,7 +34,6 @@ KRATOS_DEFINE_APPLICATION_VARIABLE( CONVECTION_DIFFUSION_APPLICATION, double, ME
 KRATOS_DEFINE_APPLICATION_VARIABLE( CONVECTION_DIFFUSION_APPLICATION, double, MEAN_VEL_OVER_ELEM_SIZE)
 KRATOS_DEFINE_APPLICATION_VARIABLE( CONVECTION_DIFFUSION_APPLICATION, double, MELT_TEMPERATURE_1)
 KRATOS_DEFINE_APPLICATION_VARIABLE( CONVECTION_DIFFUSION_APPLICATION, double, MELT_TEMPERATURE_2)
-KRATOS_DEFINE_APPLICATION_VARIABLE( CONVECTION_DIFFUSION_APPLICATION, double, PENALTY_DIRICHLET)
 KRATOS_DEFINE_APPLICATION_VARIABLE( CONVECTION_DIFFUSION_APPLICATION, double, EMBEDDED_SCALAR)
 KRATOS_DEFINE_APPLICATION_VARIABLE( CONVECTION_DIFFUSION_APPLICATION, double, PROJECTED_SCALAR1)
 KRATOS_DEFINE_APPLICATION_VARIABLE( CONVECTION_DIFFUSION_APPLICATION, double, TRANSFER_COEFFICIENT)

applications/ConvectionDiffusionApplication/custom_elements/embedded_laplacian_element.cpp

@@ -103,7 +103,7 @@ void EmbeddedLaplacianElement<TTDim>::CalculateLocalSystem(
 
         // Calculate and add local system for the positive side of the element
         AddPositiveElementSide(rLeftHandSideMatrix, rRightHandSideVector, rCurrentProcessInfo, data);
-        
+
         // Calculate and add interface terms
         AddPositiveInterfaceTerms(rLeftHandSideMatrix, rRightHandSideVector, rCurrentProcessInfo, data);
 
@@ -217,28 +217,28 @@ void EmbeddedLaplacianElement<TTDim>::AddPositiveElementSide(
         temp[n] = r_geom[n].GetSolutionStepValue(r_unknown_var);
     }
 
-    // Iterate over the positive side volume integration points 
+    // Iterate over the positive side volume integration points
     // = number of integration points * number of subdivisions on positive side of element
     const std::size_t number_of_positive_gauss_points = rData.PositiveSideWeights.size();
     for (std::size_t g = 0; g < number_of_positive_gauss_points; ++g) {
 
-        const auto& N = row(rData.PositiveSideN, g); 
+        const auto& N = row(rData.PositiveSideN, g);
         const auto& DN_DX = rData.PositiveSideDNDX[g];
-        const double weight_gauss = rData.PositiveSideWeights[g]; 
+        const double weight_gauss = rData.PositiveSideWeights[g];
 
         //Calculate the local conductivity
         const double conductivity_gauss = inner_prod(N, nodal_conductivity);
 
-        noalias(rLeftHandSideMatrix) += weight_gauss * conductivity_gauss * prod(DN_DX, trans(DN_DX)); 
+        noalias(rLeftHandSideMatrix) += weight_gauss * conductivity_gauss * prod(DN_DX, trans(DN_DX));
 
         // Calculate the local RHS (external source)
         const double q_gauss = inner_prod(N, heat_flux_local);
 
         noalias(rRightHandSideVector) += weight_gauss * q_gauss * N;
     }
-    
+
     //RHS -= K*temp
-    noalias(rRightHandSideVector) -= prod(rLeftHandSideMatrix,temp);  
+    noalias(rRightHandSideVector) -= prod(rLeftHandSideMatrix,temp);
 }
 
 template<std::size_t TTDim>
@@ -262,13 +262,13 @@ void EmbeddedLaplacianElement<TTDim>::AddPositiveInterfaceTerms(
         temp[n] = r_geom[n].GetSolutionStepValue(r_unknown_var);
     }
 
-    // Iterate over the positive side interface integration points 
+    // Iterate over the positive side interface integration points
     const std::size_t number_of_positive_gauss_points = rData.PositiveInterfaceWeights.size();
     for (std::size_t g = 0; g < number_of_positive_gauss_points; ++g) {
 
-        const auto& N = row(rData.PositiveInterfaceN, g); 
+        const auto& N = row(rData.PositiveInterfaceN, g);
         const auto& DN_DX = rData.PositiveInterfaceDNDX[g];
-        const double weight_gauss = rData.PositiveInterfaceWeights[g]; 
+        const double weight_gauss = rData.PositiveInterfaceWeights[g];
         const auto& unit_normal = rData.PositiveInterfaceUnitNormals[g];
 
         //Calculate the local conductivity
@@ -311,19 +311,19 @@ void EmbeddedLaplacianElement<TTDim>::AddNitscheBoundaryTerms(
     }
 
     // Nitsche penalty constant
-    const double gamma = rCurrentProcessInfo[PENALTY_DIRICHLET]; 
+    const double gamma = rCurrentProcessInfo[PENALTY_COEFFICIENT];
     // Dirichlet boundary value
     const double temp_bc = GetValue(EMBEDDED_SCALAR);
     // Measure of element size
     const double h = ElementSizeCalculator<TTDim,NumNodes>::AverageElementSize(r_geom);
 
-    // Iterate over the positive side interface integration points 
+    // Iterate over the positive side interface integration points
     const std::size_t number_of_positive_gauss_points = rData.PositiveInterfaceWeights.size();
     for (std::size_t g = 0; g < number_of_positive_gauss_points; ++g) {
 
-        const auto& N = row(rData.PositiveInterfaceN, g); 
+        const auto& N = row(rData.PositiveInterfaceN, g);
         const auto& DN_DX = rData.PositiveInterfaceDNDX[g];
-        const double weight_gauss = rData.PositiveInterfaceWeights[g]; 
+        const double weight_gauss = rData.PositiveInterfaceWeights[g];
         const auto& unit_normal = rData.PositiveInterfaceUnitNormals[g];
 
         //Calculate the local conductivity and temperature (at Gauss point)
@@ -341,7 +341,7 @@ void EmbeddedLaplacianElement<TTDim>::AddNitscheBoundaryTerms(
             for (std::size_t d = 0; d < TTDim; ++d) {
                 aux_bc_2 += weight_gauss * conductivity_gauss * DN_DX(i, d) * unit_normal(d);
             }
-            
+
             // Add contribution of Nitsche boundary condition to LHS
             for (std::size_t j = 0; j < NumNodes; ++j) {
                 rLeftHandSideMatrix(i, j) += aux_bc_1 * N(j) - aux_bc_2 * N(j);

applications/ConvectionDiffusionApplication/custom_python/kratos_convection_diffusion_python_application.cpp

@@ -52,7 +52,6 @@ PYBIND11_MODULE(KratosConvectionDiffusionApplication,m)
     KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, MEAN_SIZE)
     KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, MELT_TEMPERATURE_1)
     KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, MELT_TEMPERATURE_2)
-    KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, PENALTY_DIRICHLET)
     KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, EMBEDDED_SCALAR)
     KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, PROJECTED_SCALAR1)
     KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, TRANSFER_COEFFICIENT)

applications/ConvectionDiffusionApplication/tests/cpp_tests/test_embedded_laplacian_element.cpp

@@ -57,7 +57,7 @@ namespace Kratos::Testing
         const auto& r_process_info = r_test_model_part.GetProcessInfo();
 
         // Set Dirichlet penalty constant and boundary value
-        p_element->SetValue(PENALTY_DIRICHLET, 1e0);
+        p_element->SetValue(PENALTY_COEFFICIENT, 1e0);
         p_element->SetValue(EMBEDDED_SCALAR, 0.0);
 
         // Set distances for uncut element

applications/ConvectionDiffusionApplication/tests/test_convection_diffusion_embedded_solver.py

@@ -25,7 +25,7 @@ def ApplyBoundaryConditions(self):
             node.SetSolutionStepValue(KratosMultiphysics.HEAT_FLUX, 0, 1.0)
 
         # Set penalty coefficient
-        self._GetSolver().GetComputingModelPart().ProcessInfo[KratosMultiphysics.ConvectionDiffusionApplication.PENALTY_DIRICHLET] = 1.0e0
+        self._GetSolver().GetComputingModelPart().ProcessInfo[KratosMultiphysics.PENALTY_COEFFICIENT] = 1.0e0
         # Set boundary value for Nitsche imposition of DBC
         for elem in self._GetSolver().GetComputingModelPart().Elements:
             elem.SetValue(KratosMultiphysics.ConvectionDiffusionApplication.EMBEDDED_SCALAR, 0.0)

applications/DamApplication/python_scripts/dam_eigen_solver.py

@@ -1,6 +1,5 @@
 #import kratos core and applications
 import KratosMultiphysics
-import KratosMultiphysics.EigenSolversApplication as KratosEigen
 import KratosMultiphysics.StructuralMechanicsApplication as KratosStructural
 import KratosMultiphysics.DamApplication as KratosDam
 import KratosMultiphysics.PoromechanicsApplication as KratosPoro
@@ -69,7 +68,7 @@ def AddVariables(self):
 
 
     def GetMinimumBufferSize(self):
-        return 2;
+        return 2
 
 
     def AddDofs(self):

applications/OptimizationApplication/custom_utilities/response/linear_strain_energy_response_utils.cpp

@@ -156,7 +156,7 @@ void LinearStrainEnergyResponseUtils::CalculateStrainEnergyEntitySemiAnalyticSha
     Vector& rX,
     Vector& rRefRHS,
     Vector& rPerturbedRHS,
-    typename TEntityType::Pointer& pThreadLocalEntity,
+    Node::Pointer& pThreadLocalNode,
     ModelPart& rModelPart,
     const double Delta,
     const Variable<array_1d<double, 3>>& rOutputGradientVariable)
@@ -170,71 +170,59 @@ void LinearStrainEnergyResponseUtils::CalculateStrainEnergyEntitySemiAnalyticSha
 
         rEntity.GetValuesVector(rX);
 
+        rX /= 2.0;
+
         // calculate the reference value
         rEntity.CalculateRightHandSide(rRefRHS, r_process_info);
 
-        // initialize dummy element for parallelized perturbation based sensitivity calculation
-        // in each thread seperately
-        if (!pThreadLocalEntity) {
-            GeometryType::PointsArrayType nodes;
-
-            for (IndexType i = 0; i < r_geometry.size(); ++i) {
-                const auto& r_coordinates = r_geometry[i].Coordinates();
-                auto p_new_node = Kratos::make_intrusive<Node>(i+1, r_coordinates[0], r_coordinates[1], r_coordinates[2]);
-                p_new_node->SetSolutionStepVariablesList(rModelPart.pGetNodalSolutionStepVariablesList());
-                p_new_node->SetBufferSize(rModelPart.GetBufferSize());
-                nodes.push_back(p_new_node);
-            }
-
-            pThreadLocalEntity = rEntity.Create(1, nodes, rEntity.pGetProperties());
-            pThreadLocalEntity->Initialize(r_process_info);
-            pThreadLocalEntity->InitializeSolutionStep(r_process_info);
-        }
-
-        *pThreadLocalEntity = static_cast<const TEntityType&>(rEntity);
-
-        // TODO: For some reason, assignment operator of the element
-        //       assigns current position to initial position of the lhs element
-        //       hence these are corrected manually. But need to be checked
-        //       in the element/node/geometry assignment operators.
-        pThreadLocalEntity->GetGeometry().SetData(r_geometry.GetData());
-        pThreadLocalEntity->SetFlags(rEntity.GetFlags());
-        for (IndexType i = 0; i < r_geometry.size(); ++i) {
-            auto& r_node = pThreadLocalEntity->GetGeometry()[i];
-            const auto& r_orig_node = r_geometry[i];
-            r_node = r_orig_node;
+        // initialize dummy node for parallelized perturbation based sensitivity calculation
+        // in each thread separately
+        if (!pThreadLocalNode) {
+            pThreadLocalNode = Kratos::make_intrusive<Node>(1, 0, 0, 0);
         }
 
         // now calculate perturbed
         for (IndexType i = 0; i < r_geometry.size(); ++i) {
             auto& r_orig_node_sensitivity = r_geometry[i].GetValue(rOutputGradientVariable);
 
-            auto& r_node = pThreadLocalEntity->GetGeometry()[i];
-            auto& r_coordinates = r_node.Coordinates();
-            auto& r_initial_coordintes = r_node.GetInitialPosition();
+            // get the geometry node pointer
+            auto& p_node = rEntity.GetGeometry()(i);
 
-            r_initial_coordintes[0] += Delta;
+            // now copy the node data to thread local node using the operator= in Node
+            (*pThreadLocalNode) = (*p_node);
+
+            // now swap entity node with the thread local node
+            std::swap(p_node, pThreadLocalNode);
+
+            // now do the finite difference computation.
+            auto& r_coordinates = p_node->Coordinates();
+            auto& r_initial_coordinates = p_node->GetInitialPosition();
+
+            r_initial_coordinates[0] += Delta;
             r_coordinates[0] += Delta;
-            pThreadLocalEntity->CalculateRightHandSide(rPerturbedRHS, r_process_info);
-            r_initial_coordintes[0] -= Delta;
+            rEntity.CalculateRightHandSide(rPerturbedRHS, r_process_info);
+            r_initial_coordinates[0] -= Delta;
             r_coordinates[0] -= Delta;
-            AtomicAdd<double>(r_orig_node_sensitivity[0], 0.5 * inner_prod(rX, rPerturbedRHS - rRefRHS) / Delta);
+            AtomicAdd<double>(r_orig_node_sensitivity[0], inner_prod(rX, rPerturbedRHS - rRefRHS) / Delta);
 
-            r_initial_coordintes[1] += Delta;
+            r_initial_coordinates[1] += Delta;
             r_coordinates[1] += Delta;
-            pThreadLocalEntity->CalculateRightHandSide(rPerturbedRHS, r_process_info);
-            r_initial_coordintes[1] -= Delta;
+            rEntity.CalculateRightHandSide(rPerturbedRHS, r_process_info);
+            r_initial_coordinates[1] -= Delta;
             r_coordinates[1] -= Delta;
-            AtomicAdd<double>(r_orig_node_sensitivity[1], 0.5 * inner_prod(rX, rPerturbedRHS - rRefRHS) / Delta);
+            AtomicAdd<double>(r_orig_node_sensitivity[1], inner_prod(rX, rPerturbedRHS - rRefRHS) / Delta);
 
             if (domain_size == 3) {
-                r_initial_coordintes[2] += Delta;
+                r_initial_coordinates[2] += Delta;
                 r_coordinates[2] += Delta;
-                pThreadLocalEntity->CalculateRightHandSide(rPerturbedRHS, r_process_info);
-                r_initial_coordintes[2] -= Delta;
+                rEntity.CalculateRightHandSide(rPerturbedRHS, r_process_info);
+                r_initial_coordinates[2] -= Delta;
                 r_coordinates[2] -= Delta;
-                AtomicAdd<double>(r_orig_node_sensitivity[2], 0.5 * inner_prod(rX, rPerturbedRHS - rRefRHS) / Delta);
+                AtomicAdd<double>(r_orig_node_sensitivity[2], inner_prod(rX, rPerturbedRHS - rRefRHS) / Delta);
             }
+
+            // revert back the node change.
+            std::swap(p_node, pThreadLocalNode);
         }
     }
 
@@ -248,13 +236,11 @@ void LinearStrainEnergyResponseUtils::CalculateStrainEnergySemiAnalyticShapeGrad
 {
     KRATOS_TRY
 
-    using tls_element_type = std::tuple<Vector, Vector, Vector, Element::Pointer>;
-
-    using tls_condition_type = std::tuple<Vector, Vector, Vector, Condition::Pointer>;
+    using tls_type = std::tuple<Vector, Vector, Vector, Node::Pointer>;
 
     VariableUtils().SetNonHistoricalVariableToZero(rOutputGradientVariable, rModelPart.Nodes());
 
-    block_for_each(rModelPart.Elements(), tls_element_type(), [&](auto& rElement, tls_element_type& rTLS) {
+    block_for_each(rModelPart.Elements(), tls_type(), [&](auto& rElement, tls_type& rTLS) {
         CalculateStrainEnergyEntitySemiAnalyticShapeGradient(
             rElement,
             std::get<0>(rTLS),
@@ -266,7 +252,7 @@ void LinearStrainEnergyResponseUtils::CalculateStrainEnergySemiAnalyticShapeGrad
             rOutputGradientVariable);
     });
 
-    block_for_each(rModelPart.Conditions(), tls_condition_type(), [&](auto& rCondition, tls_condition_type& rTLS) {
+    block_for_each(rModelPart.Conditions(), tls_type(), [&](auto& rCondition, tls_type& rTLS) {
         CalculateStrainEnergyEntitySemiAnalyticShapeGradient(
             rCondition,
             std::get<0>(rTLS),

applications/OptimizationApplication/custom_utilities/response/linear_strain_energy_response_utils.h

@@ -76,7 +76,7 @@ class KRATOS_API(OPTIMIZATION_APPLICATION) LinearStrainEnergyResponseUtils
         Vector& rX,
         Vector& rRefRHS,
         Vector& rPerturbedRHS,
-        typename TEntityType::Pointer& pThreadLocalEntity,
+        Node::Pointer& pThreadLocalEntity,
         ModelPart& rModelPart,
         const double Delta,
         const Variable<array_1d<double, 3>>& rOutputGradientVariable);