diff --git a/applications/MeshingApplication/custom_processes/metrics_hessian_process.cpp b/applications/MeshingApplication/custom_processes/metrics_hessian_process.cpp index 2c27d7400da2..ba7d90ec3f08 100755 --- a/applications/MeshingApplication/custom_processes/metrics_hessian_process.cpp +++ b/applications/MeshingApplication/custom_processes/metrics_hessian_process.cpp @@ -20,157 +20,108 @@ namespace Kratos { -template -ComputeHessianSolMetricProcess::ComputeHessianSolMetricProcess( - ModelPart& rThisModelPart, - TVarType& rVariable, - Parameters ThisParameters - ):mThisModelPart(rThisModelPart), - mVariable(rVariable) +ComputeHessianSolMetricProcess::ComputeHessianSolMetricProcess( + ModelPart& rThisModelPart, + Variable& rVariable, + Parameters ThisParameters + ):mThisModelPart(rThisModelPart) { - Parameters default_parameters = Parameters(R"( - { - "minimal_size" : 0.1, - "maximal_size" : 10.0, - "enforce_current" : true, - "hessian_strategy_parameters": - { - "metric_variable" : ["DISTANCE"], - "estimate_interpolation_error" : false, - "interpolation_error" : 1.0e-6, - "mesh_dependent_constant" : 0.28125 - }, - "anisotropy_remeshing" : true, - "anisotropy_parameters": - { - "reference_variable_name" : "DISTANCE", - "hmin_over_hmax_anisotropic_ratio" : 1.0, - "boundary_layer_max_distance" : 1.0, - "interpolation" : "Linear" - } - })" ); + // We push the list of double variables + mrOriginVariableDoubleList.push_back(&rVariable); + // We check the parameters + Parameters default_parameters = GetDefaultParameters(); ThisParameters.RecursivelyValidateAndAssignDefaults(default_parameters); + InitializeVariables(ThisParameters); +} - mMinSize = ThisParameters["minimal_size"].GetDouble(); - mMaxSize = ThisParameters["maximal_size"].GetDouble(); - mEnforceCurrent = ThisParameters["enforce_current"].GetBool(); +/***********************************************************************************/ +/***********************************************************************************/ - // In case we have isotropic remeshing (default values) - if (ThisParameters["anisotropy_remeshing"].GetBool() == false) { - mEstimateInterpError = default_parameters["hessian_strategy_parameters"]["estimate_interpolation_error"].GetBool(); - mInterpError = default_parameters["hessian_strategy_parameters"]["interpolation_error"].GetDouble(); - mMeshConstant = default_parameters["hessian_strategy_parameters"]["mesh_dependent_constant"].GetDouble(); - mRatioReferenceVariable = default_parameters["anisotropy_parameters"]["reference_variable_name"].GetString(); - mAnisotropicRatio = default_parameters["anisotropy_parameters"]["hmin_over_hmax_anisotropic_ratio"].GetDouble(); - mBoundLayer = default_parameters["anisotropy_parameters"]["boundary_layer_max_distance"].GetDouble(); - mInterpolation = ConvertInter(default_parameters["anisotropy_parameters"]["interpolation"].GetString()); - } else { - mEstimateInterpError = ThisParameters["hessian_strategy_parameters"]["estimate_interpolation_error"].GetBool(); - mInterpError = ThisParameters["hessian_strategy_parameters"]["interpolation_error"].GetDouble(); - mMeshConstant = ThisParameters["hessian_strategy_parameters"]["mesh_dependent_constant"].GetDouble(); - mRatioReferenceVariable = ThisParameters["anisotropy_parameters"]["reference_variable_name"].GetString(); - mAnisotropicRatio = ThisParameters["anisotropy_parameters"]["hmin_over_hmax_anisotropic_ratio"].GetDouble(); - mBoundLayer = ThisParameters["anisotropy_parameters"]["boundary_layer_max_distance"].GetDouble(); - mInterpolation = ConvertInter(ThisParameters["anisotropy_parameters"]["interpolation"].GetString()); - } +ComputeHessianSolMetricProcess::ComputeHessianSolMetricProcess( + ModelPart& rThisModelPart, + ComponentType& rVariable, + Parameters ThisParameters + ):mThisModelPart(rThisModelPart) +{ + // We push the components list + mrOriginVariableComponentsList.push_back(&rVariable); + + // We check the parameters + Parameters default_parameters = GetDefaultParameters(); + ThisParameters.RecursivelyValidateAndAssignDefaults(default_parameters); + InitializeVariables(ThisParameters); } /***********************************************************************************/ /***********************************************************************************/ -template -void ComputeHessianSolMetricProcess::Execute() +void ComputeHessianSolMetricProcess::Execute() { - // Iterate in the nodes - NodesArrayType& nodes_array = mThisModelPart.Nodes(); - const int num_nodes = nodes_array.end() - nodes_array.begin(); - + // Computing auxiliar Hessian CalculateAuxiliarHessian(); // Some checks - VariableUtils().CheckVariableExists(mVariable, nodes_array); - for (auto& i_node : nodes_array) - KRATOS_ERROR_IF_NOT(i_node.Has(NODAL_H)) << "NODAL_H must be computed" << std::endl; - - // Ratio reference variable - KRATOS_ERROR_IF_NOT(KratosComponents>::Has(mRatioReferenceVariable)) << "Variable " << mRatioReferenceVariable << " is not a double variable" << std::endl; - const auto& reference_var = KratosComponents>::Get(mRatioReferenceVariable); - - // Tensor variable definition - const Variable& tensor_variable = KratosComponents>::Get("METRIC_TENSOR_"+std::to_string(TDim)+"D"); - - // Setting metric in case not defined - if (!nodes_array.begin()->Has(tensor_variable)) { - // Declaring auxiliar vector - const TensorArrayType aux_zero_vector = ZeroVector(3 * (TDim - 1)); - #pragma omp parallel for - for(int i = 0; i < num_nodes; ++i) { - auto it_node = nodes_array.begin() + i; - it_node->SetValue(tensor_variable, aux_zero_vector); - } + NodesArrayType& nodes_array = mThisModelPart.Nodes(); + if (mrOriginVariableDoubleList.size() > 0) { + VariableUtils().CheckVariableExists(*mrOriginVariableDoubleList[0], nodes_array); + } else { + VariableUtils().CheckVariableExists(*mrOriginVariableComponentsList[0], nodes_array); } - #pragma omp parallel for - for(int i = 0; i < num_nodes; ++i) { - auto it_node = nodes_array.begin() + i; - - const Vector& hessian = it_node->GetValue(AUXILIAR_HESSIAN); - - const double nodal_h = it_node->GetValue(NODAL_H); - - double element_min_size = mMinSize; - if ((element_min_size > nodal_h) && mEnforceCurrent) element_min_size = nodal_h; - double element_max_size = mMaxSize; - if ((element_max_size > nodal_h) && mEnforceCurrent) element_max_size = nodal_h; - - // Isotropic by default - double ratio = 1.0; - - if (it_node->SolutionStepsDataHas(reference_var)) { - const double ratio_reference = it_node->FastGetSolutionStepValue(reference_var); - ratio = CalculateAnisotropicRatio(ratio_reference, mAnisotropicRatio, mBoundLayer, mInterpolation); - } - - // For postprocess pourposes - it_node->SetValue(ANISOTROPIC_RATIO, ratio); - - // We compute the metric - KRATOS_DEBUG_ERROR_IF_NOT(it_node->Has(tensor_variable)) << "METRIC_TENSOR_" + std::to_string(TDim) + "D not defined for node " << it_node->Id() << std::endl; - TensorArrayType& metric = it_node->GetValue(tensor_variable); + // Checking NODAL_H + for (const auto& i_node : nodes_array) + KRATOS_ERROR_IF_NOT(i_node.Has(NODAL_H)) << "NODAL_H must be computed" << std::endl; - const double norm_metric = norm_2(metric); - if (norm_metric > 0.0) {// NOTE: This means we combine differents metrics, at the same time means that the metric should be reseted each time - const TensorArrayType& old_metric = it_node->GetValue(tensor_variable); - const TensorArrayType& new_metric = ComputeHessianMetricTensor(hessian, ratio, element_min_size, element_max_size); + // Getting dimension + const std::size_t dimension = mThisModelPart.GetProcessInfo()[DOMAIN_SIZE]; - metric = MetricsMathUtils::IntersectMetrics(old_metric, new_metric); - } else { - metric = ComputeHessianMetricTensor(hessian, ratio, element_min_size, element_max_size); - } + // Computing metric + if (dimension == 2) { // 2D + CalculateMetric<2>(); + } else if (dimension == 3) { // 3D + CalculateMetric<3>(); + } else { + KRATOS_ERROR << "Dimension can be only 2D or 3D. Dimension: " << dimension << std::endl; } } /***********************************************************************************/ /***********************************************************************************/ -template -array_1d ComputeHessianSolMetricProcess::ComputeHessianMetricTensor( +template +array_1d ComputeHessianSolMetricProcess::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 ) { + /// The type of array considered for the tensor + typedef typename std::conditional, array_1d>::type TensorArrayType; + + /// Matrix type definition + typedef BoundedMatrix MatrixType; + // We first transform the Hessian into a matrix const MatrixType hessian_matrix = MathUtils::VectorToSymmetricTensor(rHessian); // Calculating Metric parameters + const double mesh_constant = TDim == 3 ? 9.0/32.0 : 2.0/9.0; double interpolation_error = mInterpError; - if (mEstimateInterpError) - interpolation_error = 2.0/9.0 * MathUtils::Max(ElementMaxSize, ElementMaxSize * norm_frobenius(hessian_matrix)); + if (mEstimateInterpError) { + interpolation_error = mesh_constant * MathUtils::Max(ElementMaxSize, ElementMaxSize * norm_frobenius(hessian_matrix)); + } - KRATOS_ERROR_IF(interpolation_error < std::numeric_limits::epsilon()) << "ERROR: YOUR INTERPOLATION ERROR IS NEAR ZERO: " << interpolation_error << std::endl; + // We check is the interpolation error is near zero. If it is we will correct it + if (interpolation_error < std::numeric_limits::epsilon()) { + KRATOS_WARNING("ComputeHessianSolMetricProcess") << "WARNING: Your interpolation error is near zero: " << interpolation_error << ". Computing a local L(inf) upper bound of the interpolation error"<< std::endl; + const double l_square = std::pow(NodalH, 2); + for (IndexType i = 0; i < 3 * (TDim - 1); ++i) { + interpolation_error = mesh_constant * MathUtils::Max(interpolation_error, l_square * std::abs(rHessian[i])); + } + } const double c_epslilon = mMeshConstant/interpolation_error; const double min_ratio = 1.0/(ElementMinSize * ElementMinSize); const double max_ratio = 1.0/(ElementMaxSize * ElementMaxSize); @@ -219,103 +170,157 @@ array_1d ComputeHessianSolMetricProcess: /***********************************************************************************/ /***********************************************************************************/ -template -void ComputeHessianSolMetricProcess::CalculateAuxiliarHessian() +void ComputeHessianSolMetricProcess::CalculateAuxiliarHessian() { - // Iterate in the nodes - NodesArrayType& nodes_array = mThisModelPart.Nodes(); - const int num_nodes = nodes_array.end() - nodes_array.begin(); - - // Declaring auxiliar vector - const TensorArrayType aux_zero_vector = ZeroVector(3 * (TDim - 1)); - - #pragma omp parallel for - for(int i = 0; i < num_nodes; ++i) - (nodes_array.begin() + i)->SetValue(AUXILIAR_HESSIAN, aux_zero_vector); - - // Compute auxiliar gradient - ComputeNodalGradientProcess gradient_process = ComputeNodalGradientProcess(mThisModelPart, mVariable, AUXILIAR_GRADIENT, NODAL_AREA); - gradient_process.Execute(); - - // Iterate in the conditions + // Iterate in the elements ElementsArrayType& elements_array = mThisModelPart.Elements(); - const int num_elements = elements_array.end() - elements_array.begin(); + const int num_elements = static_cast(elements_array.size()); + const auto& it_element_begin = elements_array.begin(); - #pragma omp parallel for - for(int i = 0; i < num_elements; ++i) { + // Geometry information + const std::size_t dimension = mThisModelPart.GetProcessInfo()[DOMAIN_SIZE]; - Element::GeometryType& geom = (elements_array.begin() + i)->GetGeometry(); + // Declaring auxiliar vector + const Vector aux_zero_hessian = ZeroVector(3 * (dimension - 1)); + const array_1d aux_zero_vector = ZeroVector(3); - double volume; - if (geom.GetGeometryType() == GeometryData::KratosGeometryType::Kratos_Triangle2D3) { - BoundedMatrix DN_DX; - array_1d N; + // Iterate in the nodes + NodesArrayType& nodes_array = mThisModelPart.Nodes(); + const int num_nodes = static_cast(nodes_array.size()); - GeometryUtils::CalculateGeometryData(geom, DN_DX, N, volume); + // Initialize auxiliar variables + const auto& it_nodes_begin = nodes_array.begin(); + #pragma omp parallel for + for(int i_node = 0; i_node < num_nodes; ++i_node) { + auto it_node = it_nodes_begin + i_node; + it_node->SetValue(NODAL_AREA, 0.0); + it_node->SetValue(AUXILIAR_HESSIAN, aux_zero_hessian); + it_node->SetValue(AUXILIAR_GRADIENT, aux_zero_vector); + } - BoundedMatrix values; - for(IndexType i_node = 0; i_node < 3; ++i_node) { - const array_1d& aux_grad = geom[i_node].GetValue(AUXILIAR_GRADIENT); - values(i_node, 0) = aux_grad[0]; - values(i_node, 1) = aux_grad[1]; - } + // Compute auxiliar gradient + if (mrOriginVariableDoubleList.size() > 0) { + auto gradient_process = ComputeNodalGradientProcess(mThisModelPart, *mrOriginVariableDoubleList[0], AUXILIAR_GRADIENT, NODAL_AREA); + gradient_process.Execute(); + } else { + auto gradient_process = ComputeNodalGradientProcess(mThisModelPart, *mrOriginVariableComponentsList[0], AUXILIAR_GRADIENT, NODAL_AREA); + gradient_process.Execute(); + } - const BoundedMatrix& hessian = prod(trans(DN_DX), values); - const array_1d& hessian_cond = MathUtils::StressTensorToVector, array_1d>(hessian); + // Auxiliar containers + Matrix DN_DX, J0; + Vector N; + + #pragma omp parallel for firstprivate(DN_DX, N, J0) + for(int i_elem = 0; i_elem < num_elements; ++i_elem) { + auto it_elem = it_element_begin + i_elem; + auto& r_geometry = it_elem->GetGeometry(); + + // Current geometry information + const std::size_t local_space_dimension = r_geometry.LocalSpaceDimension(); + const std::size_t number_of_nodes = r_geometry.PointsNumber(); + + // Resize if needed + if (DN_DX.size1() != number_of_nodes || DN_DX.size2() != dimension) + DN_DX.resize(number_of_nodes, dimension); + if (N.size() != number_of_nodes) + N.resize(number_of_nodes); + if (J0.size1() != dimension || J0.size2() != local_space_dimension) + J0.resize(dimension, local_space_dimension); + + // The integration points + const auto& integration_method = r_geometry.GetDefaultIntegrationMethod(); + const auto& integration_points = r_geometry.IntegrationPoints(integration_method); + const std::size_t number_of_integration_points = integration_points.size(); + + // The containers of the shape functions and the local gradients + const auto& rNcontainer = r_geometry.ShapeFunctionsValues(integration_method); + const auto& rDN_DeContainer = r_geometry.ShapeFunctionsLocalGradients(integration_method); + + // 2D case + if (dimension == 2) { + for ( IndexType point_number = 0; point_number < number_of_integration_points; ++point_number ) { + // Getting the shape functions + noalias(N) = row(rNcontainer, point_number); + + // Getting the jacobians and local gradients + GeometryUtils::JacobianOnInitialConfiguration(r_geometry, integration_points[point_number], J0); + double detJ0; + Matrix InvJ0; + MathUtils::InvertMatrix(J0, InvJ0, detJ0); + const Matrix& rDN_De = rDN_DeContainer[point_number]; + GeometryUtils::ShapeFunctionsGradients(rDN_De, InvJ0, DN_DX); + + const double gauss_point_volume = integration_points[point_number].Weight() * detJ0; + + Matrix values(number_of_nodes, 2); + for(IndexType i_node = 0; i_node < number_of_nodes; ++i_node) { + const array_1d& aux_grad = r_geometry[i_node].GetValue(AUXILIAR_GRADIENT); + for (IndexType i_dim = 0; i_dim < 2; ++i_dim) + values(i_node, i_dim) = aux_grad[i_dim]; + } - for(IndexType i_node = 0; i_node < geom.size(); ++i_node) { - auto& aux_hessian = geom[i_node].GetValue(AUXILIAR_HESSIAN); - for(IndexType k = 0; k < 3; ++k) { - double& val = aux_hessian[k]; + const BoundedMatrix& hessian = prod(trans(DN_DX), values); + const array_1d& hessian_cond = MathUtils::StressTensorToVector, array_1d>(hessian); - #pragma omp atomic - val += N[i_node] * volume * hessian_cond[k]; + for(IndexType i_node = 0; i_node < number_of_nodes; ++i_node) { + auto& aux_hessian = r_geometry[i_node].GetValue(AUXILIAR_HESSIAN); + for(IndexType k = 0; k < 3; ++k) { + #pragma omp atomic + aux_hessian[k] += N[i_node] * gauss_point_volume * hessian_cond[k]; + } } } - } - else if (geom.GetGeometryType() == GeometryData::KratosGeometryType::Kratos_Tetrahedra3D4) { - BoundedMatrix DN_DX; - array_1d N; - - GeometryUtils::CalculateGeometryData(geom, DN_DX, N, volume); - - BoundedMatrix values; - for(IndexType i_node = 0; i_node < 4; ++i_node) { - const array_1d& aux_grad = geom[i_node].GetValue(AUXILIAR_GRADIENT); - values(i_node, 0) = aux_grad[0]; - values(i_node, 1) = aux_grad[1]; - values(i_node, 2) = aux_grad[2]; - } - - const BoundedMatrix hessian = prod(trans(DN_DX), values); - const array_1d& hessian_cond = MathUtils::StressTensorToVector, array_1d>(hessian); + } else { // 3D case + for ( IndexType point_number = 0; point_number < number_of_integration_points; ++point_number ) { + // Getting the shape functions + noalias(N) = row(rNcontainer, point_number); + + // Getting the jacobians and local gradients + GeometryUtils::JacobianOnInitialConfiguration(r_geometry, integration_points[point_number], J0); + double detJ0; + Matrix InvJ0; + MathUtils::InvertMatrix(J0, InvJ0, detJ0); + const Matrix& rDN_De = rDN_DeContainer[point_number]; + GeometryUtils::ShapeFunctionsGradients(rDN_De, InvJ0, DN_DX); + + const double gauss_point_volume = integration_points[point_number].Weight() * detJ0; + + Matrix values(number_of_nodes, 3); + for(IndexType i_node = 0; i_node < number_of_nodes; ++i_node) { + const array_1d& aux_grad = r_geometry[i_node].GetValue(AUXILIAR_GRADIENT); + for (IndexType i_dim = 0; i_dim < 3; ++i_dim) + values(i_node, i_dim) = aux_grad[i_dim]; + } - for(IndexType i_node = 0; i_node < geom.size(); ++i_node) { - auto& aux_hessian = geom[i_node].GetValue(AUXILIAR_HESSIAN); - for(IndexType k = 0; k < 6; ++k) { - double& val = aux_hessian[k]; + const BoundedMatrix hessian = prod(trans(DN_DX), values); + const array_1d& hessian_cond = MathUtils::StressTensorToVector, array_1d>(hessian); - #pragma omp atomic - val += N[i_node] * volume * hessian_cond[k]; + for(IndexType i_node = 0; i_node < number_of_nodes; ++i_node) { + auto& aux_hessian = r_geometry[i_node].GetValue(AUXILIAR_HESSIAN); + for(IndexType k = 0; k < 6; ++k) { + #pragma omp atomic + aux_hessian[k] += N[i_node] * gauss_point_volume * hessian_cond[k]; + } } } } - else - KRATOS_ERROR << "WARNING: YOU CAN USE JUST 2D TRIANGLES OR 3D TETRAEDRA RIGHT NOW IN THE GEOMETRY UTILS: " << geom.size() << std::endl; } #pragma omp parallel for - for(int i = 0; i < num_nodes; ++i) { - auto it_node = nodes_array.begin() + i; - it_node->GetValue(AUXILIAR_HESSIAN) /= it_node->GetValue(NODAL_AREA); + for(int i_node = 0; i_node < num_nodes; ++i_node) { + auto it_node = nodes_array.begin() + i_node; + const double nodal_area = it_node->GetValue(NODAL_AREA); + if (nodal_area > std::numeric_limits::epsilon()) { + it_node->GetValue(AUXILIAR_HESSIAN) /= nodal_area; + } } } /***********************************************************************************/ /***********************************************************************************/ -template -double ComputeHessianSolMetricProcess::CalculateAnisotropicRatio( +double ComputeHessianSolMetricProcess::CalculateAnisotropicRatio( const double Distance, const double AnisotropicRatio, const double BoundLayer, @@ -343,9 +348,144 @@ double ComputeHessianSolMetricProcess::CalculateAnisotropicRatio /***********************************************************************************/ /***********************************************************************************/ -template class ComputeHessianSolMetricProcess<2, Variable>; -template class ComputeHessianSolMetricProcess<3, Variable>; -template class ComputeHessianSolMetricProcess<2, ComponentType>; -template class ComputeHessianSolMetricProcess<3, ComponentType>; +template +void ComputeHessianSolMetricProcess::CalculateMetric() +{ + /// The type of array considered for the tensor + typedef typename std::conditional, array_1d>::type TensorArrayType; + + // Iterate in the nodes + NodesArrayType& nodes_array = mThisModelPart.Nodes(); + const int num_nodes = static_cast(nodes_array.size()); + + // Tensor variable definition + const Variable& tensor_variable = KratosComponents>::Get("METRIC_TENSOR_"+std::to_string(TDim)+"D"); + + // Setting metric in case not defined + const auto it_node_begin = nodes_array.begin(); + if (!it_node_begin->Has(tensor_variable)) { + // Declaring auxiliar vector + const TensorArrayType aux_zero_vector = ZeroVector(3 * (TDim - 1)); + VariableUtils().SetNonHistoricalVariable(tensor_variable, aux_zero_vector, nodes_array); + } + + // Ratio reference variable + KRATOS_ERROR_IF_NOT(KratosComponents>::Has(mRatioReferenceVariable)) << "Variable " << mRatioReferenceVariable << " is not a double variable" << std::endl; + const auto& reference_var = KratosComponents>::Get(mRatioReferenceVariable); + + #pragma omp parallel for + for(int i = 0; i < num_nodes; ++i) { + auto it_node = it_node_begin + i; + + const Vector& hessian = it_node->GetValue(AUXILIAR_HESSIAN); + + const double nodal_h = it_node->GetValue(NODAL_H); + + double element_min_size = mMinSize; + if ((element_min_size > nodal_h) && mEnforceCurrent) element_min_size = nodal_h; + double element_max_size = mMaxSize; + if ((element_max_size > nodal_h) && mEnforceCurrent) element_max_size = nodal_h; + + // Isotropic by default + double ratio = 1.0; + + if (it_node->SolutionStepsDataHas(reference_var)) { + const double ratio_reference = it_node->FastGetSolutionStepValue(reference_var); + ratio = CalculateAnisotropicRatio(ratio_reference, mAnisotropicRatio, mBoundLayer, mInterpolation); + } + + // For postprocess pourposes + it_node->SetValue(ANISOTROPIC_RATIO, ratio); + + // We compute the metric + KRATOS_DEBUG_ERROR_IF_NOT(it_node->Has(tensor_variable)) << "METRIC_TENSOR_" + std::to_string(TDim) + "D not defined for node " << it_node->Id() << std::endl; + TensorArrayType& metric = it_node->GetValue(tensor_variable); + + const double norm_metric = norm_2(metric); + if (norm_metric > 0.0) {// NOTE: This means we combine differents metrics, at the same time means that the metric should be reseted each time + const TensorArrayType& old_metric = it_node->GetValue(tensor_variable); + const TensorArrayType& new_metric = ComputeHessianMetricTensor(hessian, ratio, element_min_size, element_max_size, nodal_h); + + metric = MetricsMathUtils::IntersectMetrics(old_metric, new_metric); + } else { + metric = ComputeHessianMetricTensor(hessian, ratio, element_min_size, element_max_size, nodal_h); + } + } +} + +/***********************************************************************************/ +/***********************************************************************************/ + +Parameters ComputeHessianSolMetricProcess::GetDefaultParameters() const +{ + Parameters default_parameters = Parameters(R"( + { + "minimal_size" : 0.1, + "maximal_size" : 10.0, + "enforce_current" : true, + "hessian_strategy_parameters": + { + "metric_variable" : ["DISTANCE"], + "estimate_interpolation_error" : false, + "interpolation_error" : 1.0e-6, + "mesh_dependent_constant" : 0.28125 + }, + "anisotropy_remeshing" : true, + "anisotropy_parameters": + { + "reference_variable_name" : "DISTANCE", + "hmin_over_hmax_anisotropic_ratio" : 1.0, + "boundary_layer_max_distance" : 1.0, + "interpolation" : "Linear" + } + })" ); + + // Identify the dimension first + const SizeType dimension = mThisModelPart.GetProcessInfo()[DOMAIN_SIZE]; + + // The mesh dependent constant depends on dimension + if (dimension == 2) { + default_parameters["hessian_strategy_parameters"]["mesh_dependent_constant"].SetDouble(2.0/9.0); + } else if (dimension == 3) { + default_parameters["hessian_strategy_parameters"]["mesh_dependent_constant"].SetDouble(9.0/32.0); + } else { + KRATOS_ERROR << "Dimension can be only 2D or 3D. Dimension: " << dimension << std::endl; + } + + return default_parameters; +} + +/***********************************************************************************/ +/***********************************************************************************/ + +void ComputeHessianSolMetricProcess::InitializeVariables(Parameters ThisParameters) +{ + // Get default variables + Parameters default_parameters = GetDefaultParameters(); + + // Set variables + mMinSize = ThisParameters["minimal_size"].GetDouble(); + mMaxSize = ThisParameters["maximal_size"].GetDouble(); + mEnforceCurrent = ThisParameters["enforce_current"].GetBool(); + + // In case we have isotropic remeshing (default values) + if (ThisParameters["anisotropy_remeshing"].GetBool() == false) { + mEstimateInterpError = default_parameters["hessian_strategy_parameters"]["estimate_interpolation_error"].GetBool(); + mInterpError = default_parameters["hessian_strategy_parameters"]["interpolation_error"].GetDouble(); + mMeshConstant = default_parameters["hessian_strategy_parameters"]["mesh_dependent_constant"].GetDouble(); + mRatioReferenceVariable = default_parameters["anisotropy_parameters"]["reference_variable_name"].GetString(); + mAnisotropicRatio = default_parameters["anisotropy_parameters"]["hmin_over_hmax_anisotropic_ratio"].GetDouble(); + mBoundLayer = default_parameters["anisotropy_parameters"]["boundary_layer_max_distance"].GetDouble(); + mInterpolation = ConvertInter(default_parameters["anisotropy_parameters"]["interpolation"].GetString()); + } else { + mEstimateInterpError = ThisParameters["hessian_strategy_parameters"]["estimate_interpolation_error"].GetBool(); + mInterpError = ThisParameters["hessian_strategy_parameters"]["interpolation_error"].GetDouble(); + mMeshConstant = ThisParameters["hessian_strategy_parameters"]["mesh_dependent_constant"].GetDouble(); + mRatioReferenceVariable = ThisParameters["anisotropy_parameters"]["reference_variable_name"].GetString(); + mAnisotropicRatio = ThisParameters["anisotropy_parameters"]["hmin_over_hmax_anisotropic_ratio"].GetDouble(); + mBoundLayer = ThisParameters["anisotropy_parameters"]["boundary_layer_max_distance"].GetDouble(); + mInterpolation = ConvertInter(ThisParameters["anisotropy_parameters"]["interpolation"].GetString()); + } +} };// namespace Kratos. diff --git a/applications/MeshingApplication/custom_processes/metrics_hessian_process.h b/applications/MeshingApplication/custom_processes/metrics_hessian_process.h index 019a2b6fd1b8..1eb908e20cb5 100755 --- a/applications/MeshingApplication/custom_processes/metrics_hessian_process.h +++ b/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 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, array_1d>::type TensorArrayType; - - /// Matrix type definition - typedef BoundedMatrix 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& 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*> mrOriginVariableDoubleList; /// The scalar variable list to compute + std::vector 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 array_1d 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 + 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 inline std::istream& operator >> (std::istream& rIStream, - ComputeHessianSolMetricProcess& rThis); + ComputeHessianSolMetricProcess& rThis); /// output stream function -template inline std::ostream& operator << (std::ostream& rOStream, - const ComputeHessianSolMetricProcess& rThis) + const ComputeHessianSolMetricProcess& rThis) { rThis.PrintInfo(rOStream); rOStream << std::endl; diff --git a/applications/MeshingApplication/custom_python/add_processes_to_python.cpp b/applications/MeshingApplication/custom_python/add_processes_to_python.cpp index d17d5c2dc593..2f82a9f4b33f 100644 --- a/applications/MeshingApplication/custom_python/add_processes_to_python.cpp +++ b/applications/MeshingApplication/custom_python/add_processes_to_python.cpp @@ -89,27 +89,18 @@ void AddProcessesToPython(pybind11::module& m) .def(py::init>, Parameters>()) ; - // HESSIAN DOUBLE - py::class_>, ComputeHessianSolMetricProcess<2, Variable>::Pointer, Process>(m, "ComputeHessianSolMetricProcess2D") + // HESSIAN PROCESS + py::class_(m, "ComputeHessianSolMetricProcess") .def(py::init&>()) .def(py::init&, Parameters>()) - ; - - py::class_>, ComputeHessianSolMetricProcess<3, Variable>::Pointer, Process>(m, "ComputeHessianSolMetricProcess3D") - .def(py::init&>()) - .def(py::init&, Parameters>()) - ; - - // HESSIAN ARRAY 1D - py::class_, ComputeHessianSolMetricProcess<2, ComponentType>::Pointer, Process>(m, "ComputeHessianSolMetricProcessComp2D") .def(py::init()) .def(py::init()) ; - py::class_, ComputeHessianSolMetricProcess<3, ComponentType>::Pointer, Process>(m, "ComputeHessianSolMetricProcessComp3D") - .def(py::init()) - .def(py::init()) - ; + 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>::Pointer, Process>(m, "MetricErrorProcess2D") diff --git a/applications/MeshingApplication/python_scripts/mmg_process.py b/applications/MeshingApplication/python_scripts/mmg_process.py index 5b6ad2b34a8f..3b938d33d270 100644 --- a/applications/MeshingApplication/python_scripts/mmg_process.py +++ b/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() diff --git a/applications/MeshingApplication/tests/cpp_tests/test_metric_process.cpp b/applications/MeshingApplication/tests/cpp_tests/test_metric_process.cpp index 6a7886146fdc..ad813121427f 100644 --- a/applications/MeshingApplication/tests/cpp_tests/test_metric_process.cpp +++ b/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, Historical> GradientType; + typedef ComputeNodalGradientProcess 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, Historical> GradientType; + typedef ComputeNodalGradientProcess 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> hessian_process = ComputeHessianSolMetricProcess<2, Variable>(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> hessian_process = ComputeHessianSolMetricProcess<3, Variable>(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); diff --git a/applications/MeshingApplication/tests/mmg_lagrangian_test/beam2D_line_load_test_parameters.json b/applications/MeshingApplication/tests/mmg_lagrangian_test/beam2D_line_load_test_parameters.json index e2a763f92225..2458860d740d 100644 --- a/applications/MeshingApplication/tests/mmg_lagrangian_test/beam2D_line_load_test_parameters.json +++ b/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", diff --git a/applications/MeshingApplication/tests/mmg_lagrangian_test/beam2D_line_load_test_result.sol b/applications/MeshingApplication/tests/mmg_lagrangian_test/beam2D_line_load_test_result.sol index 52f10242c6e4..11bede126df0 100644 --- a/applications/MeshingApplication/tests/mmg_lagrangian_test/beam2D_line_load_test_result.sol +++ b/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 diff --git a/applications/MeshingApplication/tests/mmg_lagrangian_test/beam2D_test_parameters.json b/applications/MeshingApplication/tests/mmg_lagrangian_test/beam2D_test_parameters.json index 187c3394ffd0..430c1f67a432 100644 --- a/applications/MeshingApplication/tests/mmg_lagrangian_test/beam2D_test_parameters.json +++ b/applications/MeshingApplication/tests/mmg_lagrangian_test/beam2D_test_parameters.json @@ -86,16 +86,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", diff --git a/applications/MeshingApplication/tests/mmg_lagrangian_test/beam2D_test_result.sol b/applications/MeshingApplication/tests/mmg_lagrangian_test/beam2D_test_result.sol index ff352c6eec42..e6e4dc52bc20 100644 --- a/applications/MeshingApplication/tests/mmg_lagrangian_test/beam2D_test_result.sol +++ b/applications/MeshingApplication/tests/mmg_lagrangian_test/beam2D_test_result.sol @@ -7,77 +7,77 @@ Dimension SolAtVertices 205 1 3 - 35.9101637379947 0 35.9101637379947 - 37.4367319907306 0 37.4367319907306 - 26.2566573854866 0 26.2566573854866 - 33.410918127346 0 33.410918127346 + 28.3767908497022 0 28.3767908497022 + 29.5822299755478 0 29.5822299755478 + 20.7475871590569 0 20.7475871590569 + 26.400746180205 0 26.400746180205 30.2668616137121 0 30.2668616137121 - 39.4995524908847 0 39.4995524908847 - 21.7518979323714 0 21.7518979323714 + 31.2128905550287 0 31.2128905550287 + 19.4301897481549 0 19.4301897481549 30.2668616137121 0 30.2668616137121 - 34.7899307544646 0 34.7899307544646 - 37.2336861625557 0 37.2336861625557 - 36.0263728933714 0 36.0263728933714 - 25.0209084013989 0 25.0209084013989 - 33.2528253574337 0 33.2528253574337 - 32.7410519353801 0 32.7410519353801 - 32.3652989726995 0 32.3652989726995 - 34.6480898721677 0 34.6480898721677 - 30.4481389105709 0 30.4481389105709 - 35.7736766879115 0 35.7736766879115 - 34.5222528912967 0 34.5222528912967 - 38.5869366910055 0 38.5869366910055 - 31.8229556412735 0 31.8229556412735 - 28.8360160894753 0 28.8360160894753 - 32.8305343431724 0 32.8305343431724 - 35.2779086996321 0 35.2779086996321 - 27.5604426825636 0 27.5604426825636 - 31.369834994686 0 31.369834994686 - 29.7131102638918 0 29.7131102638918 - 27.2120608961529 0 27.2120608961529 - 32.3129342013113 0 32.3129342013113 - 27.0409428791352 0 27.0409428791352 - 28.3068155213527 0 28.3068155213527 - 27.879208923983 0 27.879208923983 - 29.6339719883714 0 29.6339719883714 - 26.3374102967325 0 26.3374102967325 - 25.084751519167 0 25.084751519167 - 27.7256981050374 0 27.7256981050374 - 25.9932252275203 0 25.9932252275203 - 25.5242330886323 0 25.5242330886323 - 24.3569152890673 0 24.3569152890673 - 23.5078082385583 0 23.5078082385583 - 25.5003560190692 0 25.5003560190692 - 24.6357056743481 0 24.6357056743481 - 24.2516297234798 0 24.2516297234798 - 23.2508413731498 0 23.2508413731498 - 23.9496149419695 0 23.9496149419695 - 22.695867515866 0 22.695867515866 - 22.868593567692 0 22.868593567692 - 22.4404691219244 0 22.4404691219244 - 21.8135888475102 0 21.8135888475102 - 22.2702038034889 0 22.2702038034889 - 21.3863489488884 0 21.3863489488884 - 21.3263847639816 0 21.3263847639816 - 20.7023246411292 0 20.7023246411292 - 19.9829308632535 0 19.9829308632535 - 20.6919954413751 0 20.6919954413751 - 20.1100772583546 0 20.1100772583546 - 19.3895747360983 0 19.3895747360983 - 19.3479987996405 0 19.3479987996405 - 19.1049270364835 0 19.1049270364835 - 19.2993646488359 0 19.2993646488359 - 18.7602306557437 0 18.7602306557437 - 19.4308765674892 0 19.4308765674892 - 17.8599794183255 0 17.8599794183255 - 17.9198645601064 0 17.9198645601064 - 17.9608749777344 0 17.9608749777344 - 17.4075287793985 0 17.4075287793985 - 17.8079007695172 0 17.8079007695172 - 16.657995600629 0 16.657995600629 - 16.6834949119661 0 16.6834949119661 - 16.9050455076346 0 16.9050455076346 - 16.5444845098659 0 16.5444845098659 + 27.4893622951255 0 27.4893622951255 + 29.4191756831891 0 29.4191756831891 + 28.4691132075417 0 28.4691132075417 + 19.7681024622638 0 19.7681024622638 + 26.2717884449097 0 26.2717884449097 + 25.8676674158485 0 25.8676674158485 + 25.5717623624025 0 25.5717623624025 + 27.3731367364782 0 27.3731367364782 + 24.0614418225759 0 24.0614418225759 + 28.2644400405265 0 28.2644400405265 + 27.2744115348091 0 27.2744115348091 + 30.4848745867609 0 30.4848745867609 + 25.1417487980258 0 25.1417487980258 + 22.7828130998142 0 22.7828130998142 + 25.9408670860144 0 25.9408670860144 + 27.8735919013541 0 27.8735919013541 + 21.7754650238298 0 21.7754650238298 + 24.7871275245345 0 24.7871275245345 + 23.4793981846999 0 23.4793981846999 + 21.5011709128867 0 21.5011709128867 + 25.5313929672512 0 25.5313929672512 + 21.3646835331044 0 21.3646835331044 + 22.3672846498199 0 22.3672846498199 + 22.0305998650072 0 22.0305998650072 + 23.4105083749785 0 23.4105083749785 + 20.8080417793664 0 20.8080417793664 + 19.8237995614721 0 19.8237995614721 + 21.9027498606528 0 21.9027498606528 + 20.5395133343945 0 20.5395133343945 + 20.169959338492 0 20.169959338492 + 19.2430320996938 0 19.2430320996938 + 18.5777491249299 0 18.5777491249299 + 20.1445280898505 0 20.1445280898505 + 19.4671189628613 0 19.4671189628613 + 19.1645956066766 0 19.1645956066766 + 18.3745189015379 0 18.3745189015379 + 18.9192505950814 0 18.9192505950814 + 17.9362305271928 0 17.9362305271928 + 18.0710245414658 0 18.0710245414658 + 17.7336008635128 0 17.7336008635128 + 17.238901676324 0 17.238901676324 + 17.5923621873972 0 17.5923621873972 + 16.9014981887464 0 16.9014981887464 + 16.8525941384334 0 16.8525941384334 + 16.3544100295369 0 16.3544100295369 + 16.2448121764785 0 16.2448121764785 + 16.3510343200608 0 16.3510343200608 + 16.2448121764785 0 16.2448121764785 + 16.2448121764785 0 16.2448121764785 + 16.2448121764785 0 16.2448121764785 + 16.2448121764785 0 16.2448121764785 + 16.2448121764785 0 16.2448121764785 + 16.2448121764785 0 16.2448121764785 + 16.2448121764785 0 16.2448121764785 + 16.2448121764785 0 16.2448121764785 + 16.2448121764785 0 16.2448121764785 + 16.2448121764785 0 16.2448121764785 + 16.2448121764785 0 16.2448121764785 + 16.2448121764785 0 16.2448121764785 + 16.2448121764785 0 16.2448121764785 + 16.2448121764785 0 16.2448121764785 + 16.2448121764785 0 16.2448121764785 + 16.2448121764785 0 16.2448121764785 16.2448121764785 0 16.2448121764785 16.2448121764785 0 16.2448121764785 16.2448121764785 0 16.2448121764785 diff --git a/applications/MeshingApplication/tests/mmg_lagrangian_test/beam3D_test_result.sol b/applications/MeshingApplication/tests/mmg_lagrangian_test/beam3D_test_result.sol index 493d92a0bd65..31e3e0c64943 100644 --- a/applications/MeshingApplication/tests/mmg_lagrangian_test/beam3D_test_result.sol +++ b/applications/MeshingApplication/tests/mmg_lagrangian_test/beam3D_test_result.sol @@ -7,138 +7,138 @@ Dimension SolAtVertices 224 1 3 -27.7035626161233 0 27.7035626161233 0 0 27.7035626161233 -17.022367691514 0 17.022367691514 0 0 17.022367691514 -25.355481664933 0 25.355481664933 0 0 25.355481664933 -20.0288569979345 0 20.0288569979345 0 0 20.0288569979345 -21.89484350378 0 21.89484350378 0 0 21.89484350378 -22.5371603355171 0 22.5371603355171 0 0 22.5371603355171 -25.2328942120125 0 25.2328942120125 0 0 25.2328942120125 -24.6906878741862 0 24.6906878741862 0 0 24.6906878741862 -23.5209851368877 0 23.5209851368877 0 0 23.5209851368877 -27.2051885688751 0 27.2051885688751 0 0 27.2051885688751 -23.7038914066742 0 23.7038914066742 0 0 23.7038914066742 -26.8232916744992 0 26.8232916744992 0 0 26.8232916744992 -20.5420941066609 0 20.5420941066609 0 0 20.5420941066609 -20.9751300360015 0 20.9751300360015 0 0 20.9751300360015 -30.0155679986682 0 30.0155679986682 0 0 30.0155679986682 -23.4678846419606 0 23.4678846419606 0 0 23.4678846419606 -15.2378861117701 0 15.2378861117701 0 0 15.2378861117701 -24.4602537914722 0 24.4602537914722 0 0 24.4602537914722 -26.0312536439935 0 26.0312536439935 0 0 26.0312536439935 -26.3754649289435 0 26.3754649289435 0 0 26.3754649289435 -25.3588617914023 0 25.3588617914023 0 0 25.3588617914023 -25.8562007727243 0 25.8562007727243 0 0 25.8562007727243 -25.9820866445989 0 25.9820866445989 0 0 25.9820866445989 -28.0485555161631 0 28.0485555161631 0 0 28.0485555161631 -28.2067151374578 0 28.2067151374578 0 0 28.2067151374578 -28.0943518256723 0 28.0943518256723 0 0 28.0943518256723 -24.6232520628818 0 24.6232520628818 0 0 24.6232520628818 -23.1627548939625 0 23.1627548939625 0 0 23.1627548939625 -25.0746783214843 0 25.0746783214843 0 0 25.0746783214843 -27.0514757048651 0 27.0514757048651 0 0 27.0514757048651 -28.5184631116839 0 28.5184631116839 0 0 28.5184631116839 -21.5808089749154 0 21.5808089749154 0 0 21.5808089749154 -22.3410583905441 0 22.3410583905441 0 0 22.3410583905441 -22.0333439029875 0 22.0333439029875 0 0 22.0333439029875 -21.7450316639833 0 21.7450316639833 0 0 21.7450316639833 -23.1482758220858 0 23.1482758220858 0 0 23.1482758220858 -20.5911061371322 0 20.5911061371322 0 0 20.5911061371322 -19.3067696885139 0 19.3067696885139 0 0 19.3067696885139 -20.0891252144247 0 20.0891252144247 0 0 20.0891252144247 -21.5935804930098 0 21.5935804930098 0 0 21.5935804930098 -21.3806124773869 0 21.3806124773869 0 0 21.3806124773869 -22.9581844181729 0 22.9581844181729 0 0 22.9581844181729 -20.8562612449711 0 20.8562612449711 0 0 20.8562612449711 -19.5486594086639 0 19.5486594086639 0 0 19.5486594086639 -18.2028063118509 0 18.2028063118509 0 0 18.2028063118509 -20.2619604911674 0 20.2619604911674 0 0 20.2619604911674 -18.7699614525379 0 18.7699614525379 0 0 18.7699614525379 -18.2346814216115 0 18.2346814216115 0 0 18.2346814216115 -19.2139607699384 0 19.2139607699384 0 0 19.2139607699384 -16.6312785781068 0 16.6312785781068 0 0 16.6312785781068 -20.7226895955545 0 20.7226895955545 0 0 20.7226895955545 -19.7766515181522 0 19.7766515181522 0 0 19.7766515181522 -17.1565673048547 0 17.1565673048547 0 0 17.1565673048547 -17.6439367464756 0 17.6439367464756 0 0 17.6439367464756 -16.1116279129533 0 16.1116279129533 0 0 16.1116279129533 -19.2843584261176 0 19.2843584261176 0 0 19.2843584261176 -16.967474634563 0 16.967474634563 0 0 16.967474634563 -16.9902676810897 0 16.9902676810897 0 0 16.9902676810897 -17.5921558562483 0 17.5921558562483 0 0 17.5921558562483 -14.8276810794386 0 14.8276810794386 0 0 14.8276810794386 -18.2855785011589 0 18.2855785011589 0 0 18.2855785011589 -14.8661717707709 0 14.8661717707709 0 0 14.8661717707709 -14.7686109193122 0 14.7686109193122 0 0 14.7686109193122 -16.9554609692555 0 16.9554609692555 0 0 16.9554609692555 -13.455540017933 0 13.455540017933 0 0 13.455540017933 -14.185483630148 0 14.185483630148 0 0 14.185483630148 -15.8277064200399 0 15.8277064200399 0 0 15.8277064200399 -14.1849599388393 0 14.1849599388393 0 0 14.1849599388393 -14.1018061112717 0 14.1018061112717 0 0 14.1018061112717 -13.4442013135134 0 13.4442013135134 0 0 13.4442013135134 -13.5999714770016 0 13.5999714770016 0 0 13.5999714770016 -14.5481936084951 0 14.5481936084951 0 0 14.5481936084951 -14.5479848494158 0 14.5479848494158 0 0 14.5479848494158 -12.9909304345219 0 12.9909304345219 0 0 12.9909304345219 -12.8092665161218 0 12.8092665161218 0 0 12.8092665161218 -11.3983751171588 0 11.3983751171588 0 0 11.3983751171588 -11.6144472344102 0 11.6144472344102 0 0 11.6144472344102 -12.9055162943599 0 12.9055162943599 0 0 12.9055162943599 -11.9333785334854 0 11.9333785334854 0 0 11.9333785334854 -12.1137609258182 0 12.1137609258182 0 0 12.1137609258182 -12.5696566494541 0 12.5696566494541 0 0 12.5696566494541 -11.3550189726419 0 11.3550189726419 0 0 11.3550189726419 -10.3677942410406 0 10.3677942410406 0 0 10.3677942410406 -10.1398620246949 0 10.1398620246949 0 0 10.1398620246949 -10.760157630112 0 10.760157630112 0 0 10.760157630112 -10.259438150323 0 10.259438150323 0 0 10.259438150323 -9.65921823604312 0 9.65921823604312 0 0 9.65921823604312 -9.1418657201665 0 9.1418657201665 0 0 9.1418657201665 -9.39233794545191 0 9.39233794545191 0 0 9.39233794545191 -9.02694252618019 0 9.02694252618019 0 0 9.02694252618019 -9.72916342998384 0 9.72916342998384 0 0 9.72916342998384 -9.58652539038531 0 9.58652539038531 0 0 9.58652539038531 -8.32678915384327 0 8.32678915384327 0 0 8.32678915384327 -9.27454215520464 0 9.27454215520464 0 0 9.27454215520464 -9.53277137698498 0 9.53277137698498 0 0 9.53277137698498 -8.89894818733828 0 8.89894818733828 0 0 8.89894818733828 -8.51660823683054 0 8.51660823683054 0 0 8.51660823683054 -8.22609130240348 0 8.22609130240348 0 0 8.22609130240348 -7.87642627781283 0 7.87642627781283 0 0 7.87642627781283 -8.33949942560789 0 8.33949942560789 0 0 8.33949942560789 -8.5807679085591 0 8.5807679085591 0 0 8.5807679085591 -7.85107300473627 0 7.85107300473627 0 0 7.85107300473627 -8.42153377613511 0 8.42153377613511 0 0 8.42153377613511 -6.91498232488896 0 6.91498232488896 0 0 6.91498232488896 -8.02491969678893 0 8.02491969678893 0 0 8.02491969678893 -7.40582739573172 0 7.40582739573172 0 0 7.40582739573172 -7.15130340264934 0 7.15130340264934 0 0 7.15130340264934 +27.7057416888572 0 27.7057416888572 0 0 27.7057416888572 +17.0223617137705 0 17.0223617137705 0 0 17.0223617137705 +25.3568602712687 0 25.3568602712687 0 0 25.3568602712687 +20.0283499660681 0 20.0283499660681 0 0 20.0283499660681 +21.8959336459052 0 21.8959336459052 0 0 21.8959336459052 +22.5376651387608 0 22.5376651387608 0 0 22.5376651387608 +25.2348849289845 0 25.2348849289845 0 0 25.2348849289845 +24.6906579607744 0 24.6906579607744 0 0 24.6906579607744 +23.521905668729 0 23.521905668729 0 0 23.521905668729 +27.2028032365598 0 27.2028032365598 0 0 27.2028032365598 +23.7054976030952 0 23.7054976030952 0 0 23.7054976030952 +26.8243396718601 0 26.8243396718601 0 0 26.8243396718601 +20.5413014388128 0 20.5413014388128 0 0 20.5413014388128 +20.9738129236526 0 20.9738129236526 0 0 20.9738129236526 +30.0139579495844 0 30.0139579495844 0 0 30.0139579495844 +23.4691868016015 0 23.4691868016015 0 0 23.4691868016015 +15.237420818953 0 15.237420818953 0 0 15.237420818953 +24.4591378882655 0 24.4591378882655 0 0 24.4591378882655 +26.0316962124953 0 26.0316962124953 0 0 26.0316962124953 +26.3766673195715 0 26.3766673195715 0 0 26.3766673195715 +25.3604866056377 0 25.3604866056377 0 0 25.3604866056377 +25.8577215644443 0 25.8577215644443 0 0 25.8577215644443 +25.981061888996 0 25.981061888996 0 0 25.981061888996 +28.0462801370087 0 28.0462801370087 0 0 28.0462801370087 +28.2077696756407 0 28.2077696756407 0 0 28.2077696756407 +28.0932587116364 0 28.0932587116364 0 0 28.0932587116364 +24.621674048973 0 24.621674048973 0 0 24.621674048973 +23.1612958766903 0 23.1612958766903 0 0 23.1612958766903 +25.076709398064 0 25.076709398064 0 0 25.076709398064 +27.0524684691216 0 27.0524684691216 0 0 27.0524684691216 +28.5193560954211 0 28.5193560954211 0 0 28.5193560954211 +21.5825402240725 0 21.5825402240725 0 0 21.5825402240725 +22.3433381659784 0 22.3433381659784 0 0 22.3433381659784 +22.0344961916228 0 22.0344961916228 0 0 22.0344961916228 +21.7435600014087 0 21.7435600014087 0 0 21.7435600014087 +23.1466536412384 0 23.1466536412384 0 0 23.1466536412384 +20.5931898199178 0 20.5931898199178 0 0 20.5931898199178 +19.3092821821558 0 19.3092821821558 0 0 19.3092821821558 +20.0869812521228 0 20.0869812521228 0 0 20.0869812521228 +21.5948417598932 0 21.5948417598932 0 0 21.5948417598932 +21.3821103133146 0 21.3821103133146 0 0 21.3821103133146 +22.9607448032594 0 22.9607448032594 0 0 22.9607448032594 +20.8536951624291 0 20.8536951624291 0 0 20.8536951624291 +19.5487537117075 0 19.5487537117075 0 0 19.5487537117075 +18.2031009644258 0 18.2031009644258 0 0 18.2031009644258 +20.2648693637971 0 20.2648693637971 0 0 20.2648693637971 +18.7672124285481 0 18.7672124285481 0 0 18.7672124285481 +18.2374163066341 0 18.2374163066341 0 0 18.2374163066341 +19.2112877416363 0 19.2112877416363 0 0 19.2112877416363 +16.6338551874949 0 16.6338551874949 0 0 16.6338551874949 +20.7195954174141 0 20.7195954174141 0 0 20.7195954174141 +19.7739110349597 0 19.7739110349597 0 0 19.7739110349597 +17.154178028332 0 17.154178028332 0 0 17.154178028332 +17.6463152926316 0 17.6463152926316 0 0 17.6463152926316 +16.1089946721111 0 16.1089946721111 0 0 16.1089946721111 +19.2868107746792 0 19.2868107746792 0 0 19.2868107746792 +16.9649846850946 0 16.9649846850946 0 0 16.9649846850946 +16.9927382086824 0 16.9927382086824 0 0 16.9927382086824 +17.5944569502259 0 17.5944569502259 0 0 17.5944569502259 +14.8299139455351 0 14.8299139455351 0 0 14.8299139455351 +18.2829290757951 0 18.2829290757951 0 0 18.2829290757951 +14.8636697773977 0 14.8636697773977 0 0 14.8636697773977 +14.7707289214118 0 14.7707289214118 0 0 14.7707289214118 +16.9581109949642 0 16.9581109949642 0 0 16.9581109949642 +13.4534439547323 0 13.4534439547323 0 0 13.4534439547323 +14.1832885653418 0 14.1832885653418 0 0 14.1832885653418 +15.8298718306675 0 15.8298718306675 0 0 15.8298718306675 +14.1868526744314 0 14.1868526744314 0 0 14.1868526744314 +14.1039576683544 0 14.1039576683544 0 0 14.1039576683544 +13.4464138623935 0 13.4464138623935 0 0 13.4464138623935 +13.59766283544 0 13.59766283544 0 0 13.59766283544 +14.5457976841306 0 14.5457976841306 0 0 14.5457976841306 +14.5454286035634 0 14.5454286035634 0 0 14.5454286035634 +12.9888387133305 0 12.9888387133305 0 0 12.9888387133305 +12.8115209520957 0 12.8115209520957 0 0 12.8115209520957 +11.4004674578534 0 11.4004674578534 0 0 11.4004674578534 +11.6164851666254 0 11.6164851666254 0 0 11.6164851666254 +12.9031208486348 0 12.9031208486348 0 0 12.9031208486348 +11.9356126768213 0 11.9356126768213 0 0 11.9356126768213 +12.1115039293692 0 12.1115039293692 0 0 12.1115039293692 +12.5672877919797 0 12.5672877919797 0 0 12.5672877919797 +11.3571498489514 0 11.3571498489514 0 0 11.3571498489514 +10.3658906660789 0 10.3658906660789 0 0 10.3658906660789 +10.1377788249695 0 10.1377788249695 0 0 10.1377788249695 +10.7620991175534 0 10.7620991175534 0 0 10.7620991175534 +10.2613172442114 0 10.2613172442114 0 0 10.2613172442114 +9.65746213268743 0 9.65746213268743 0 0 9.65746213268743 +9.14200957455898 0 9.14200957455898 0 0 9.14200957455898 +9.39404070277063 0 9.39404070277063 0 0 9.39404070277063 +9.02493721337671 0 9.02493721337671 0 0 9.02493721337671 +9.73107079638935 0 9.73107079638935 0 0 9.73107079638935 +9.58463971047051 0 9.58463971047051 0 0 9.58463971047051 +8.32532694867183 0 8.32532694867183 0 0 8.32532694867183 +9.27274720466372 0 9.27274720466372 0 0 9.27274720466372 +9.53087104007251 0 9.53087104007251 0 0 9.53087104007251 +8.89728994493967 0 8.89728994493967 0 0 8.89728994493967 +8.51814956619836 0 8.51814956619836 0 0 8.51814956619836 +8.2277167139307 0 8.2277167139307 0 0 8.2277167139307 +7.87801847229764 0 7.87801847229764 0 0 7.87801847229764 +8.3377890067277 0 8.3377890067277 0 0 8.3377890067277 +8.57897370257185 0 8.57897370257185 0 0 8.57897370257185 +7.85272438427157 0 7.85272438427157 0 0 7.85272438427157 +8.41978512232926 0 8.41978512232926 0 0 8.41978512232926 +6.91637815456035 0 6.91637815456035 0 0 6.91637815456035 +8.02629418003215 0 8.02629418003215 0 0 8.02629418003215 +7.40436757626627 0 7.40436757626627 0 0 7.40436757626627 +7.15278252449708 0 7.15278252449708 0 0 7.15278252449708 7.21509701078706 0 7.21509701078706 0 0 7.21509701078706 -6.96124373392133 0 6.96124373392133 0 0 6.96124373392133 -6.77481781258337 0 6.77481781258337 0 0 6.77481781258337 -6.41416375479511 0 6.41416375479511 0 0 6.41416375479511 -7.4135495366291 0 7.4135495366291 0 0 7.4135495366291 -6.40460849965883 0 6.40460849965883 0 0 6.40460849965883 +6.96249910773807 0 6.96249910773807 0 0 6.96249910773807 +6.77365040868535 0 6.77365040868535 0 0 6.77365040868535 +6.41549949953902 0 6.41549949953902 0 0 6.41549949953902 +7.41354953662909 0 7.41354953662909 0 0 7.41354953662909 +6.40341568349489 0 6.40341568349489 0 0 6.40341568349489 7.21509701078706 0 7.21509701078706 0 0 7.21509701078706 -6.0844104143319 0 6.0844104143319 0 0 6.0844104143319 -8.81187620865784 0 8.81187620865784 0 0 8.81187620865784 -7.4135495366291 0 7.4135495366291 0 0 7.4135495366291 -6.67848542327411 0 6.67848542327411 0 0 6.67848542327411 +6.08441041433191 0 6.08441041433191 0 0 6.08441041433191 +8.81187620865783 0 8.81187620865783 0 0 8.81187620865783 +7.41354953662909 0 7.41354953662909 0 0 7.41354953662909 +6.67848542327412 0 6.67848542327412 0 0 6.67848542327412 5.53505395873405 0 5.53505395873405 0 0 5.53505395873405 -9.3422758815462 0 9.3422758815462 0 0 9.3422758815462 -7.53158400810561 0 7.53158400810561 0 0 7.53158400810561 +9.34227588154621 0 9.34227588154621 0 0 9.34227588154621 +7.5315840081056 0 7.5315840081056 0 0 7.5315840081056 8.82872300048485 0 8.82872300048485 0 0 8.82872300048485 7.2232577301389 0 7.2232577301389 0 0 7.2232577301389 -9.3422758815462 0 9.3422758815462 0 0 9.3422758815462 +9.34227588154621 0 9.34227588154621 0 0 9.34227588154621 8.82872300048485 0 8.82872300048485 0 0 8.82872300048485 5.53505395873405 0 5.53505395873405 0 0 5.53505395873405 -6.0844285066826 0 6.0844285066826 0 0 6.0844285066826 -6.59887122782258 0 6.59887122782258 0 0 6.59887122782258 +6.08442850668259 0 6.08442850668259 0 0 6.08442850668259 +6.59887122782256 0 6.59887122782256 0 0 6.59887122782256 5.53505395873405 0 5.53505395873405 0 0 5.53505395873405 6.67849296286621 0 6.67849296286621 0 0 6.67849296286621 -7.22330841065305 0 7.22330841065305 0 0 7.22330841065305 -7.22330841065305 0 7.22330841065305 0 0 7.22330841065305 +7.22330841065306 0 7.22330841065306 0 0 7.22330841065306 +7.22330841065306 0 7.22330841065306 0 0 7.22330841065306 5.53505395873405 0 5.53505395873405 0 0 5.53505395873405 5.64528686853841 0 5.64528686853841 0 0 5.64528686853841 5.53505395873405 0 5.53505395873405 0 0 5.53505395873405 @@ -146,34 +146,34 @@ SolAtVertices 5.53505395873405 0 5.53505395873405 0 0 5.53505395873405 5.53505395873405 0 5.53505395873405 0 0 5.53505395873405 5.53505395873405 0 5.53505395873405 0 0 5.53505395873405 -6.6054973151869 0 6.6054973151869 0 0 6.6054973151869 +6.60549731518692 0 6.60549731518692 0 0 6.60549731518692 5.64528686853841 0 5.64528686853841 0 0 5.64528686853841 -6.67848266745382 0 6.67848266745382 0 0 6.67848266745382 +6.67848266745381 0 6.67848266745381 0 0 6.67848266745381 5.53505395873405 0 5.53505395873405 0 0 5.53505395873405 -6.67848970352474 0 6.67848970352474 0 0 6.67848970352474 +6.67848970352473 0 6.67848970352473 0 0 6.67848970352473 6.92865943359195 0 6.92865943359195 0 0 6.92865943359195 5.53505395873405 0 5.53505395873405 0 0 5.53505395873405 7.1109948746102 0 7.1109948746102 0 0 7.1109948746102 5.53505395873405 0 5.53505395873405 0 0 5.53505395873405 7.1109948746102 0 7.1109948746102 0 0 7.1109948746102 8.77818372511676 0 8.77818372511676 0 0 8.77818372511676 -7.28428182558793 0 7.28428182558793 0 0 7.28428182558793 +7.28428182558795 0 7.28428182558795 0 0 7.28428182558795 5.54924766849245 0 5.54924766849245 0 0 5.54924766849245 -6.67848970352474 0 6.67848970352474 0 0 6.67848970352474 +6.67848970352473 0 6.67848970352473 0 0 6.67848970352473 9.52636417462794 0 9.52636417462794 0 0 9.52636417462794 5.53505395873405 0 5.53505395873405 0 0 5.53505395873405 -9.99760121892725 0 9.99760121892725 0 0 9.99760121892725 -8.88188464672214 0 8.88188464672214 0 0 8.88188464672214 -9.99760121892725 0 9.99760121892725 0 0 9.99760121892725 +9.99760121892726 0 9.99760121892726 0 0 9.99760121892726 +8.88188464672215 0 8.88188464672215 0 0 8.88188464672215 +9.99760121892726 0 9.99760121892726 0 0 9.99760121892726 9.52636417462794 0 9.52636417462794 0 0 9.52636417462794 5.53505395873405 0 5.53505395873405 0 0 5.53505395873405 -8.88188464672214 0 8.88188464672214 0 0 8.88188464672214 -6.26500647448228 0 6.26500647448228 0 0 6.26500647448228 +8.88188464672215 0 8.88188464672215 0 0 8.88188464672215 +6.26500647448227 0 6.26500647448227 0 0 6.26500647448227 6.36214889882061 0 6.36214889882061 0 0 6.36214889882061 5.53505395873405 0 5.53505395873405 0 0 5.53505395873405 6.67848696393902 0 6.67848696393902 0 0 6.67848696393902 -7.22328702528642 0 7.22328702528642 0 0 7.22328702528642 -7.22328702528642 0 7.22328702528642 0 0 7.22328702528642 +7.22328702528641 0 7.22328702528641 0 0 7.22328702528641 +7.22328702528641 0 7.22328702528641 0 0 7.22328702528641 5.53505395873405 0 5.53505395873405 0 0 5.53505395873405 5.64528459952134 0 5.64528459952134 0 0 5.64528459952134 5.53505395873405 0 5.53505395873405 0 0 5.53505395873405 @@ -188,9 +188,9 @@ SolAtVertices 5.53505395873405 0 5.53505395873405 0 0 5.53505395873405 5.53505395873405 0 5.53505395873405 0 0 5.53505395873405 5.53505395873405 0 5.53505395873405 0 0 5.53505395873405 -7.66265550765843 0 7.66265550765843 0 0 7.66265550765843 +7.66265550765842 0 7.66265550765842 0 0 7.66265550765842 5.53505395873405 0 5.53505395873405 0 0 5.53505395873405 -7.66265550765843 0 7.66265550765843 0 0 7.66265550765843 +7.66265550765842 0 7.66265550765842 0 0 7.66265550765842 6.67848759034891 0 6.67848759034891 0 0 6.67848759034891 5.53505395873405 0 5.53505395873405 0 0 5.53505395873405 5.53505395873405 0 5.53505395873405 0 0 5.53505395873405 @@ -207,15 +207,15 @@ SolAtVertices 5.61685517643149 0 5.61685517643149 0 0 5.61685517643149 11.3215991787032 0 11.3215991787032 0 0 11.3215991787032 11.3215991787032 0 11.3215991787032 0 0 11.3215991787032 -5.64971940375446 0 5.64971940375446 0 0 5.64971940375446 -5.64971940375446 0 5.64971940375446 0 0 5.64971940375446 +5.64971940375447 0 5.64971940375447 0 0 5.64971940375447 +5.64971940375447 0 5.64971940375447 0 0 5.64971940375447 5.53505395873405 0 5.53505395873405 0 0 5.53505395873405 5.53505395873405 0 5.53505395873405 0 0 5.53505395873405 5.53505395873405 0 5.53505395873405 0 0 5.53505395873405 5.81843413180583 0 5.81843413180583 0 0 5.81843413180583 5.56752500708012 0 5.56752500708012 0 0 5.56752500708012 5.53505395873405 0 5.53505395873405 0 0 5.53505395873405 -6.20162000898058 0 6.20162000898058 0 0 6.20162000898058 +6.20162000898057 0 6.20162000898057 0 0 6.20162000898057 5.6922563024583 0 5.6922563024583 0 0 5.6922563024583 5.66174406534089 0 5.66174406534089 0 0 5.66174406534089 5.53505395873405 0 5.53505395873405 0 0 5.53505395873405 diff --git a/kratos/processes/compute_nodal_gradient_process.cpp b/kratos/processes/compute_nodal_gradient_process.cpp index 6ab097fa28f6..f7f9cce07dbf 100644 --- a/kratos/processes/compute_nodal_gradient_process.cpp +++ b/kratos/processes/compute_nodal_gradient_process.cpp @@ -7,15 +7,12 @@ // License: BSD License // Kratos default license: kratos/license.txt // -// Main authors: Riccardo Rossi +// Main authors: Riccardo Rossi // Vicente Mataix Ferrandiz // // /* System includes */ -#include -#include -#include /* External includes */ @@ -26,122 +23,119 @@ namespace Kratos { -template< int TDim, class TVarType, HistoricalValues THist> -void ComputeNodalGradientProcess::Execute() +template +void ComputeNodalGradientProcess::Execute() { KRATOS_TRY; - + // Set to zero ClearGradient(); - - BoundedMatrix DN_DX; - array_1d N; - double Volume; - - #pragma omp parallel for private(DN_DX, N, Volume) - for(int i=0; i(mrModelPart.Elements().size()); ++i) { - auto it_elem = mrModelPart.ElementsBegin()+i; - Element::GeometryType& geom = it_elem->GetGeometry(); - GeometryUtils::CalculateGeometryData(geom, DN_DX, N, Volume); - - array_1d values; - for(std::size_t i=0; i grad = prod(trans(DN_DX), values); - - for(std::size_t i=0; i(mrModelPart.Elements().size()); ++i_elem) { + auto it_elem = it_element_begin + i_elem; + auto& r_geometry = it_elem->GetGeometry(); + + // Current geometry information + const std::size_t local_space_dimension = r_geometry.LocalSpaceDimension(); + const std::size_t number_of_nodes = r_geometry.PointsNumber(); + + // Resize if needed + if (DN_DX.size1() != number_of_nodes || DN_DX.size2() != dimension) + DN_DX.resize(number_of_nodes, dimension); + if (N.size() != number_of_nodes) + N.resize(number_of_nodes); + if (J0.size1() != dimension || J0.size2() != local_space_dimension) + J0.resize(dimension, local_space_dimension); + + // The integration points + const auto& integration_method = r_geometry.GetDefaultIntegrationMethod(); + const auto& integration_points = r_geometry.IntegrationPoints(integration_method); + const std::size_t number_of_integration_points = integration_points.size(); + + Vector values(number_of_nodes); + if (mrOriginVariableDoubleList.size() > 0) { + for(std::size_t i_node=0; i_node::InvertMatrix(J0, InvJ0, detJ0); + const Matrix& rDN_De = rDN_DeContainer[point_number]; + GeometryUtils::ShapeFunctionsGradients(rDN_De, InvJ0, DN_DX); + + const Vector grad = prod(trans(DN_DX), values); + const double gauss_point_volume = integration_points[point_number].Weight() * detJ0; + + for(std::size_t i_node=0; i_node& r_gradient = GetGradient(r_geometry, i_node); + for(std::size_t k=0; k -ComputeNodalGradientProcess<2, Variable, Historical>::ComputeNodalGradientProcess( - ModelPart& rModelPart, - Variable& rOriginVariable, - Variable >& rGradientVariable, - Variable& rAreaVariable) - :mrModelPart(rModelPart), mrOriginVariable(rOriginVariable), mrGradientVariable(rGradientVariable), mrAreaVariable(rAreaVariable) -{ - KRATOS_TRY - - VariableUtils().CheckVariableExists(rOriginVariable, mrModelPart.Nodes()); - VariableUtils().CheckVariableExists(rGradientVariable, mrModelPart.Nodes()); - // In case the area variable is not initialized we initialize it - auto& r_nodes = rModelPart.Nodes(); - if (!r_nodes.begin()->Has( rAreaVariable )) { - VariableUtils().SetNonHistoricalVariable(rAreaVariable, 0.0, r_nodes); - } - - KRATOS_CATCH("") -} /***********************************************************************************/ /***********************************************************************************/ template<> -ComputeNodalGradientProcess<2, Variable, NonHistorical>::ComputeNodalGradientProcess( - ModelPart& rModelPart, - Variable& rOriginVariable, - Variable >& rGradientVariable, +ComputeNodalGradientProcess::ComputeNodalGradientProcess( + ModelPart& rModelPart, + Variable& rOriginVariable, + Variable >& rGradientVariable, Variable& rAreaVariable) - :mrModelPart(rModelPart), mrOriginVariable(rOriginVariable), mrGradientVariable(rGradientVariable), mrAreaVariable(rAreaVariable) + :mrModelPart(rModelPart), mrGradientVariable(rGradientVariable), mrAreaVariable(rAreaVariable) { KRATOS_TRY - - VariableUtils().CheckVariableExists(rOriginVariable, mrModelPart.Nodes()); - // In case the area or gradient variable is not initialized we initialize it - auto& r_nodes = rModelPart.Nodes(); - if (!r_nodes.begin()->Has( rGradientVariable )) { - const array_1d zero_vector = ZeroVector(3); - VariableUtils().SetNonHistoricalVariable(rGradientVariable, zero_vector, r_nodes); - } - if (!r_nodes.begin()->Has( rAreaVariable )) { - VariableUtils().SetNonHistoricalVariable(rAreaVariable, 0.0, r_nodes); - } - - KRATOS_CATCH("") -} -/***********************************************************************************/ -/***********************************************************************************/ + // We push the list of double variables + mrOriginVariableDoubleList.push_back(&rOriginVariable); -template<> -ComputeNodalGradientProcess<3, Variable, Historical>::ComputeNodalGradientProcess( - ModelPart& rModelPart, - Variable& rOriginVariable, - Variable >& rGradientVariable, - Variable& rAreaVariable) - :mrModelPart(rModelPart), mrOriginVariable(rOriginVariable), mrGradientVariable(rGradientVariable), mrAreaVariable(rAreaVariable) -{ - KRATOS_TRY - VariableUtils().CheckVariableExists(rOriginVariable, mrModelPart.Nodes()); VariableUtils().CheckVariableExists(rGradientVariable, mrModelPart.Nodes()); - // In case the area variable is not initialized we initialize it + // In case the area or gradient variable is not initialized we initialize it auto& r_nodes = rModelPart.Nodes(); if (!r_nodes.begin()->Has( rAreaVariable )) { VariableUtils().SetNonHistoricalVariable(rAreaVariable, 0.0, r_nodes); } - + KRATOS_CATCH("") } @@ -149,15 +143,18 @@ ComputeNodalGradientProcess<3, Variable, Historical>::ComputeNodalGradie /***********************************************************************************/ template<> -ComputeNodalGradientProcess<3, Variable, NonHistorical>::ComputeNodalGradientProcess( - ModelPart& rModelPart, - Variable& rOriginVariable, - Variable >& rGradientVariable, +ComputeNodalGradientProcess::ComputeNodalGradientProcess( + ModelPart& rModelPart, + Variable& rOriginVariable, + Variable >& rGradientVariable, Variable& rAreaVariable) - :mrModelPart(rModelPart), mrOriginVariable(rOriginVariable), mrGradientVariable(rGradientVariable), mrAreaVariable(rAreaVariable) + :mrModelPart(rModelPart), mrGradientVariable(rGradientVariable), mrAreaVariable(rAreaVariable) { KRATOS_TRY - + + // We push the list of double variables + mrOriginVariableDoubleList.push_back(&rOriginVariable); + VariableUtils().CheckVariableExists(rOriginVariable, mrModelPart.Nodes()); // In case the area or gradient variable is not initialized we initialize it auto& r_nodes = rModelPart.Nodes(); @@ -168,7 +165,7 @@ ComputeNodalGradientProcess<3, Variable, NonHistorical>::ComputeNodalGra if (!r_nodes.begin()->Has( rAreaVariable )) { VariableUtils().SetNonHistoricalVariable(rAreaVariable, 0.0, r_nodes); } - + KRATOS_CATCH("") } @@ -176,50 +173,26 @@ ComputeNodalGradientProcess<3, Variable, NonHistorical>::ComputeNodalGra /***********************************************************************************/ template<> -ComputeNodalGradientProcess<2, component_type, Historical>::ComputeNodalGradientProcess( - ModelPart& rModelPart, - component_type& rOriginVariable, - Variable >& rGradientVariable, +ComputeNodalGradientProcess::ComputeNodalGradientProcess( + ModelPart& rModelPart, + ComponentType& rOriginVariable, + Variable >& rGradientVariable, Variable& rAreaVariable) - :mrModelPart(rModelPart), mrOriginVariable(rOriginVariable), mrGradientVariable(rGradientVariable), mrAreaVariable(rAreaVariable) + :mrModelPart(rModelPart), mrGradientVariable(rGradientVariable), mrAreaVariable(rAreaVariable) { KRATOS_TRY - - VariableUtils().CheckVariableExists(rOriginVariable, mrModelPart.Nodes()); - VariableUtils().CheckVariableExists(rGradientVariable, mrModelPart.Nodes()); - // In case the area variable is not initialized we initialize it - auto& r_nodes = rModelPart.Nodes(); - if (!r_nodes.begin()->Has( rAreaVariable )) { - VariableUtils().SetNonHistoricalVariable(rAreaVariable, 0.0, r_nodes); - } - - KRATOS_CATCH("") -} -/***********************************************************************************/ -/***********************************************************************************/ + // We push the components list + mrOriginVariableComponentsList.push_back(&rOriginVariable); -template<> -ComputeNodalGradientProcess<2, component_type, NonHistorical>::ComputeNodalGradientProcess( - ModelPart& rModelPart, - component_type& rOriginVariable, - Variable >& rGradientVariable, - Variable& rAreaVariable) - :mrModelPart(rModelPart), mrOriginVariable(rOriginVariable), mrGradientVariable(rGradientVariable), mrAreaVariable(rAreaVariable) -{ - KRATOS_TRY - VariableUtils().CheckVariableExists(rOriginVariable, mrModelPart.Nodes()); + VariableUtils().CheckVariableExists(rGradientVariable, mrModelPart.Nodes()); // In case the area or gradient variable is not initialized we initialize it auto& r_nodes = rModelPart.Nodes(); - if (!r_nodes.begin()->Has( rGradientVariable )) { - const array_1d zero_vector = ZeroVector(3); - VariableUtils().SetNonHistoricalVariable(rGradientVariable, zero_vector, r_nodes); - } if (!r_nodes.begin()->Has( rAreaVariable )) { VariableUtils().SetNonHistoricalVariable(rAreaVariable, 0.0, r_nodes); } - + KRATOS_CATCH("") } @@ -227,39 +200,18 @@ ComputeNodalGradientProcess<2, component_type, NonHistorical>::ComputeNodalGradi /***********************************************************************************/ template<> -ComputeNodalGradientProcess<3, component_type, Historical>::ComputeNodalGradientProcess( - ModelPart& rModelPart, - component_type& rOriginVariable, - Variable >& rGradientVariable, +ComputeNodalGradientProcess::ComputeNodalGradientProcess( + ModelPart& rModelPart, + ComponentType& rOriginVariable, + Variable >& rGradientVariable, Variable& rAreaVariable) - :mrModelPart(rModelPart), mrOriginVariable(rOriginVariable), mrGradientVariable(rGradientVariable), mrAreaVariable(rAreaVariable) + :mrModelPart(rModelPart), mrGradientVariable(rGradientVariable), mrAreaVariable(rAreaVariable) { KRATOS_TRY - - VariableUtils().CheckVariableExists(rOriginVariable, mrModelPart.Nodes()); - VariableUtils().CheckVariableExists(rGradientVariable, mrModelPart.Nodes()); - // In case the area variable is not initialized we initialize it - auto& r_nodes = rModelPart.Nodes(); - if (!r_nodes.begin()->Has( rAreaVariable )) { - VariableUtils().SetNonHistoricalVariable(rAreaVariable, 0.0, r_nodes); - } - - KRATOS_CATCH("") -} -/***********************************************************************************/ -/***********************************************************************************/ + // We push the components list + mrOriginVariableComponentsList.push_back(&rOriginVariable); -template<> -ComputeNodalGradientProcess<3, component_type, NonHistorical>::ComputeNodalGradientProcess( - ModelPart& rModelPart, - component_type& rOriginVariable, - Variable >& rGradientVariable, - Variable& rAreaVariable) - :mrModelPart(rModelPart), mrOriginVariable(rOriginVariable), mrGradientVariable(rGradientVariable), mrAreaVariable(rAreaVariable) -{ - KRATOS_TRY - VariableUtils().CheckVariableExists(rOriginVariable, mrModelPart.Nodes()); // In case the area or gradient variable is not initialized we initialize it auto& r_nodes = rModelPart.Nodes(); @@ -270,54 +222,15 @@ ComputeNodalGradientProcess<3, component_type, NonHistorical>::ComputeNodalGradi if (!r_nodes.begin()->Has( rAreaVariable )) { VariableUtils().SetNonHistoricalVariable(rAreaVariable, 0.0, r_nodes); } - - KRATOS_CATCH("") -} - -/***********************************************************************************/ -/***********************************************************************************/ - -template<> -void ComputeNodalGradientProcess<2, Variable, Historical>::ClearGradient() -{ - #pragma omp parallel for - for(int i = 0; i < static_cast(mrModelPart.Nodes().size()); ++i) { - auto it_node=mrModelPart.NodesBegin()+i; - it_node->SetValue(mrAreaVariable, 0.0); - it_node->FastGetSolutionStepValue(mrGradientVariable).clear(); - } -} - -/***********************************************************************************/ -/***********************************************************************************/ - -template<> -void ComputeNodalGradientProcess<3, Variable, Historical>::ClearGradient() -{ - #pragma omp parallel for - for(int i = 0; i < static_cast(mrModelPart.Nodes().size()); ++i) { - auto it_node=mrModelPart.NodesBegin()+i; - it_node->SetValue(mrAreaVariable, 0.0); - it_node->FastGetSolutionStepValue(mrGradientVariable).clear(); - } -} -template<> -void ComputeNodalGradientProcess<2, component_type, Historical>::ClearGradient() -{ - #pragma omp parallel for - for(int i = 0; i < static_cast(mrModelPart.Nodes().size()); ++i) { - auto it_node=mrModelPart.NodesBegin()+i; - it_node->SetValue(mrAreaVariable, 0.0); - it_node->FastGetSolutionStepValue(mrGradientVariable).clear(); - } + KRATOS_CATCH("") } /***********************************************************************************/ /***********************************************************************************/ template<> -void ComputeNodalGradientProcess<3, component_type, Historical>::ClearGradient() +void ComputeNodalGradientProcess::ClearGradient() { #pragma omp parallel for for(int i = 0; i < static_cast(mrModelPart.Nodes().size()); ++i) { @@ -331,47 +244,15 @@ void ComputeNodalGradientProcess<3, component_type, Historical>::ClearGradient() /***********************************************************************************/ template <> -void ComputeNodalGradientProcess<2, Variable, NonHistorical>::ClearGradient() -{ - const array_1d AuxZeroVector = ZeroVector(3); - - #pragma omp parallel for - for(int i = 0; i < static_cast(mrModelPart.Nodes().size()); ++i) { - auto it_node=mrModelPart.NodesBegin()+i; - it_node->SetValue(mrAreaVariable, 0.0); - it_node->SetValue(mrGradientVariable, AuxZeroVector); - } -} - -/***********************************************************************************/ -/***********************************************************************************/ - -template <> -void ComputeNodalGradientProcess<3, Variable, NonHistorical>::ClearGradient() +void ComputeNodalGradientProcess::ClearGradient() { - const array_1d AuxZeroVector = ZeroVector(3); - - #pragma omp parallel for - for(int i = 0; i < static_cast(mrModelPart.Nodes().size()); ++i) { - auto it_node=mrModelPart.NodesBegin()+i; - it_node->SetValue(mrAreaVariable, 0.0); - it_node->SetValue(mrGradientVariable, AuxZeroVector); - } -} + const array_1d aux_zero_vector = ZeroVector(3); -/***********************************************************************************/ -/***********************************************************************************/ - -template <> -void ComputeNodalGradientProcess<2, component_type, NonHistorical>::ClearGradient() -{ - const array_1d AuxZeroVector = ZeroVector(3); - #pragma omp parallel for for(int i = 0; i < static_cast(mrModelPart.Nodes().size()); ++i) { auto it_node=mrModelPart.NodesBegin()+i; it_node->SetValue(mrAreaVariable, 0.0); - it_node->SetValue(mrGradientVariable, AuxZeroVector); + it_node->SetValue(mrGradientVariable, aux_zero_vector); } } @@ -379,120 +260,12 @@ void ComputeNodalGradientProcess<2, component_type, NonHistorical>::ClearGradien /***********************************************************************************/ template <> -void ComputeNodalGradientProcess<3, component_type, NonHistorical>::ClearGradient() -{ - const array_1d AuxZeroVector = ZeroVector(3); - - #pragma omp parallel for - for(int i = 0; i < static_cast(mrModelPart.Nodes().size()); ++i) { - auto it_node=mrModelPart.NodesBegin()+i; - it_node->SetValue(mrAreaVariable, 0.0); - it_node->SetValue(mrGradientVariable, AuxZeroVector); - } -} - -/***********************************************************************************/ -/***********************************************************************************/ - -template <> -double& ComputeNodalGradientProcess<2, Variable, Historical>::GetGradient( - Element::GeometryType& rThisGeometry, - unsigned int i, - unsigned int k - ) -{ - double& val = rThisGeometry[i].FastGetSolutionStepValue(mrGradientVariable)[k]; - - return val; -} - -/***********************************************************************************/ -/***********************************************************************************/ - -template <> -double& ComputeNodalGradientProcess<3, Variable, Historical>::GetGradient( - Element::GeometryType& rThisGeometry, - unsigned int i, - unsigned int k - ) -{ - double& val = rThisGeometry[i].FastGetSolutionStepValue(mrGradientVariable)[k]; - - return val; -} - -/***********************************************************************************/ -/***********************************************************************************/ - -template <> -double& ComputeNodalGradientProcess<2, component_type, Historical>::GetGradient( - Element::GeometryType& rThisGeometry, - unsigned int i, - unsigned int k - ) -{ - double& val = rThisGeometry[i].FastGetSolutionStepValue(mrGradientVariable)[k]; - - return val; -} - -/***********************************************************************************/ -/***********************************************************************************/ - -template <> -double& ComputeNodalGradientProcess<3, component_type, Historical>::GetGradient( - Element::GeometryType& rThisGeometry, - unsigned int i, - unsigned int k - ) -{ - double& val = rThisGeometry[i].FastGetSolutionStepValue(mrGradientVariable)[k]; - - return val; -} - -/***********************************************************************************/ -/***********************************************************************************/ - -template <> -double& ComputeNodalGradientProcess<2, Variable, NonHistorical>::GetGradient( - Element::GeometryType& rThisGeometry, - unsigned int i, - unsigned int k - ) -{ - double& val = rThisGeometry[i].GetValue(mrGradientVariable)[k]; - - return val; -} - -/***********************************************************************************/ -/***********************************************************************************/ - -template <> -double& ComputeNodalGradientProcess<3, Variable, NonHistorical>::GetGradient( - Element::GeometryType& rThisGeometry, - unsigned int i, - unsigned int k - ) -{ - double& val = rThisGeometry[i].GetValue(mrGradientVariable)[k]; - - return val; -} - -/***********************************************************************************/ -/***********************************************************************************/ - -template <> -double& ComputeNodalGradientProcess<2, component_type, NonHistorical>::GetGradient( +array_1d& ComputeNodalGradientProcess::GetGradient( Element::GeometryType& rThisGeometry, - unsigned int i, - unsigned int k + unsigned int i ) { - double& val = rThisGeometry[i].GetValue(mrGradientVariable)[k]; - + array_1d& val = rThisGeometry[i].FastGetSolutionStepValue(mrGradientVariable); return val; } @@ -500,14 +273,12 @@ double& ComputeNodalGradientProcess<2, component_type, NonHistorical>::GetGradie /***********************************************************************************/ template <> -double& ComputeNodalGradientProcess<3, component_type, NonHistorical>::GetGradient( +array_1d& ComputeNodalGradientProcess::GetGradient( Element::GeometryType& rThisGeometry, - unsigned int i, - unsigned int k + unsigned int i ) { - double& val = rThisGeometry[i].GetValue(mrGradientVariable)[k]; - + array_1d& val = rThisGeometry[i].GetValue(mrGradientVariable); return val; } @@ -515,20 +286,7 @@ double& ComputeNodalGradientProcess<3, component_type, NonHistorical>::GetGradie /***********************************************************************************/ template <> -void ComputeNodalGradientProcess<2, Variable, Historical>::PonderateGradient() -{ - #pragma omp parallel for - for(int i = 0; i < static_cast(mrModelPart.Nodes().size()); ++i) { - auto it_node = mrModelPart.NodesBegin()+i; - it_node->FastGetSolutionStepValue(mrGradientVariable) /= it_node->GetValue(mrAreaVariable); - } -} - -/***********************************************************************************/ -/***********************************************************************************/ - -template <> -void ComputeNodalGradientProcess<3, Variable, Historical>::PonderateGradient() +void ComputeNodalGradientProcess::PonderateGradient() { #pragma omp parallel for for(int i = 0; i < static_cast(mrModelPart.Nodes().size()); ++i) { @@ -541,35 +299,7 @@ void ComputeNodalGradientProcess<3, Variable, Historical>::PonderateGrad /***********************************************************************************/ template <> -void ComputeNodalGradientProcess<2, component_type, Historical>::PonderateGradient() -{ - #pragma omp parallel for - for(int i = 0; i < static_cast(mrModelPart.Nodes().size()); ++i) - { - auto it_node=mrModelPart.NodesBegin()+i; - it_node->FastGetSolutionStepValue(mrGradientVariable) /= it_node->GetValue(mrAreaVariable); - } -} - -/***********************************************************************************/ -/***********************************************************************************/ - -template <> -void ComputeNodalGradientProcess<3, component_type, Historical>::PonderateGradient() -{ - #pragma omp parallel for - for(int i = 0; i < static_cast(mrModelPart.Nodes().size()); ++i) - { - auto it_node=mrModelPart.NodesBegin()+i; - it_node->FastGetSolutionStepValue(mrGradientVariable) /= it_node->GetValue(mrAreaVariable); - } -} - -/***********************************************************************************/ -/***********************************************************************************/ - -template <> -void ComputeNodalGradientProcess<2, Variable, NonHistorical>::PonderateGradient() +void ComputeNodalGradientProcess::PonderateGradient() { #pragma omp parallel for for(int i = 0; i < static_cast(mrModelPart.Nodes().size()); ++i) @@ -578,56 +308,10 @@ void ComputeNodalGradientProcess<2, Variable, NonHistorical>::PonderateG it_node->GetValue(mrGradientVariable) /= it_node->GetValue(mrAreaVariable); } } - -/***********************************************************************************/ -/***********************************************************************************/ - -template <> -void ComputeNodalGradientProcess<3, Variable, NonHistorical>::PonderateGradient() -{ - #pragma omp parallel for - for(int i = 0; i < static_cast(mrModelPart.Nodes().size()); ++i) { - auto it_node=mrModelPart.NodesBegin()+i; - it_node->GetValue(mrGradientVariable) /= it_node->GetValue(mrAreaVariable); - } -} - -/***********************************************************************************/ -/***********************************************************************************/ - -template <> -void ComputeNodalGradientProcess<2, component_type, NonHistorical>::PonderateGradient() -{ - #pragma omp parallel for - for(int i = 0; i < static_cast(mrModelPart.Nodes().size()); ++i) { - auto it_node=mrModelPart.NodesBegin()+i; - it_node->GetValue(mrGradientVariable) /= it_node->GetValue(mrAreaVariable); - } -} - -/***********************************************************************************/ -/***********************************************************************************/ - -template <> -void ComputeNodalGradientProcess<3, component_type, NonHistorical>::PonderateGradient() -{ - #pragma omp parallel for - for(int i = 0; i < static_cast(mrModelPart.Nodes().size()); ++i) { - auto it_node=mrModelPart.NodesBegin()+i; - it_node->GetValue(mrGradientVariable) /= it_node->GetValue(mrAreaVariable); - } -} - /***********************************************************************************/ /***********************************************************************************/ -template class ComputeNodalGradientProcess<2, Variable, Historical>; -template class ComputeNodalGradientProcess<2, Variable, NonHistorical>; -template class ComputeNodalGradientProcess<3, Variable, Historical>; -template class ComputeNodalGradientProcess<3, Variable, NonHistorical>; -template class ComputeNodalGradientProcess<2, component_type, Historical>; -template class ComputeNodalGradientProcess<2, component_type, NonHistorical>; -template class ComputeNodalGradientProcess<3, component_type, Historical>; -template class ComputeNodalGradientProcess<3, component_type, NonHistorical>; +template class ComputeNodalGradientProcess; +template class ComputeNodalGradientProcess; } /* namespace Kratos.*/ diff --git a/kratos/processes/compute_nodal_gradient_process.h b/kratos/processes/compute_nodal_gradient_process.h index d0f966f1f645..70b621f1045a 100644 --- a/kratos/processes/compute_nodal_gradient_process.h +++ b/kratos/processes/compute_nodal_gradient_process.h @@ -20,7 +20,6 @@ // Project includes #include "includes/define.h" -#include "includes/enums.h" #include "processes/process.h" #include "includes/model_part.h" @@ -34,12 +33,12 @@ namespace Kratos ///@name Type Definitions ///@{ - typedef VariableComponent< VectorComponentAdaptor > > component_type; + typedef VariableComponent< VectorComponentAdaptor > > ComponentType; ///@} ///@name Enum's ///@{ - + ///@} ///@name Functions ///@{ @@ -48,6 +47,16 @@ namespace Kratos ///@name Kratos Classes ///@{ +/** + * @brief This struct is used in order to identify when using the hitorical and non historical variables + */ +struct ComputeNodalGradientProcessSettings +{ + // Defining clearer options + constexpr static bool SaveAsHistoricalVariable = true; + constexpr static bool SaveAsNonHistoricalVariable = false; +}; + /** * @class ComputeNodalGradientProcess * @ingroup KratosCore @@ -55,11 +64,9 @@ namespace Kratos * @details This process computes the gradient of a certain variable stored in the nodes * @author Riccardo Rossi * @author Vicente Mataix Ferrandiz - * @tparam TDim The dimension of the problem - * @tparam TVarType The variable type - * @tparam THist If the variable is historical or not + * @tparam THistorical If the variable is historical or not */ -template< int TDim, class TVarType, HistoricalValues THist> +template class KRATOS_API(KRATOS_CORE) ComputeNodalGradientProcess : public Process { @@ -74,18 +81,21 @@ class KRATOS_API(KRATOS_CORE) ComputeNodalGradientProcess ///@name Life Cycle ///@{ - /// Default constructor. + /// Default constructor. (double) ComputeNodalGradientProcess( ModelPart& rModelPart, - TVarType& rOriginVariable, + Variable& rOriginVariable, Variable >& rGradientVariable, Variable& rAreaVariable = NODAL_AREA - ) :mrModelPart(rModelPart), - mrOriginVariable(rOriginVariable), - mrGradientVariable(rGradientVariable), - mrAreaVariable(rAreaVariable) - { - } + ); + + /// Default constructor. (component) + ComputeNodalGradientProcess( + ModelPart& rModelPart, + ComponentType& rOriginVariable, + Variable >& rGradientVariable, + Variable& rAreaVariable = NODAL_AREA + ); /// Destructor. ~ComputeNodalGradientProcess() override @@ -198,11 +208,12 @@ class KRATOS_API(KRATOS_CORE) ComputeNodalGradientProcess ///@} ///@name Member Variables ///@{ - - ModelPart& mrModelPart; // The main model part - TVarType& mrOriginVariable; // The scalar variable to compute - Variable >& mrGradientVariable; // The resultant gradient variable - Variable& mrAreaVariable; // The auxiliar area variable + + ModelPart& mrModelPart; /// The main model part + std::vector*> mrOriginVariableDoubleList; /// The scalar variable list to compute + std::vector mrOriginVariableComponentsList; /// The scalar variable list to compute (components) + Variable >& mrGradientVariable; /// The resultant gradient variable + Variable& mrAreaVariable; /// The auxiliar area variable ///@} ///@name Private Operators @@ -213,7 +224,7 @@ class KRATOS_API(KRATOS_CORE) ComputeNodalGradientProcess ///@{ // TODO: Try to use enable_if!!! - + /** * This clears the gradient */ @@ -223,14 +234,12 @@ class KRATOS_API(KRATOS_CORE) ComputeNodalGradientProcess * This gets the gradient value * @param rThisGeometry The geometry of the element * @param i The node index - * @param k The component index */ - double& GetGradient( + array_1d& GetGradient( Element::GeometryType& rThisGeometry, - unsigned int i, - unsigned int k + unsigned int i ); - + /** * This divides the gradient value by the nodal area */ @@ -274,7 +283,7 @@ class KRATOS_API(KRATOS_CORE) ComputeNodalGradientProcess /// input stream function // inline std::istream& operator >> (std::istream& rIStream, // ComputeNodalGradientProcess& rThis); -// +// // /// output stream function // inline std::ostream& operator << (std::ostream& rOStream, // const ComputeNodalGradientProcess& rThis) @@ -282,13 +291,13 @@ class KRATOS_API(KRATOS_CORE) ComputeNodalGradientProcess // rThis.PrintInfo(rOStream); // rOStream << std::endl; // rThis.PrintData(rOStream); -// +// // return rOStream; // } ///@} - + } // namespace Kratos. -#endif // KRATOS_COMPUTE_GRADIENT_PROCESS_INCLUDED defined +#endif // KRATOS_COMPUTE_GRADIENT_PROCESS_INCLUDED defined diff --git a/kratos/python/add_processes_to_python.cpp b/kratos/python/add_processes_to_python.cpp index 241f65204e21..a0896bcb9426 100644 --- a/kratos/python/add_processes_to_python.cpp +++ b/kratos/python/add_processes_to_python.cpp @@ -280,42 +280,26 @@ void AddProcessesToPython(pybind11::module& m) ; /* Historical */ - // DOUBLE - py::class_, Historical>, ComputeNodalGradientProcess<2, Variable, Historical>::Pointer, Process>(m,"ComputeNodalGradientProcess2D") - .def(py::init&, Variable >& , Variable& >()) - ; - - py::class_, Historical>, ComputeNodalGradientProcess<3, Variable, Historical>::Pointer, Process>(m,"ComputeNodalGradientProcess3D") - .def(py::init&, Variable >& , Variable& >()) - ; - - // COMPONENT - py::class_, ComputeNodalGradientProcess<2, component_type, Historical>::Pointer, Process>(m,"ComputeNodalGradientProcessComp2D") - .def(py::init >& , Variable& >()) - ; - - py::class_, ComputeNodalGradientProcess<3, component_type, Historical>::Pointer, Process>(m,"ComputeNodalGradientProcessComp3D") - .def(py::init >& , Variable& >()) - ; - + py::class_, ComputeNodalGradientProcess::Pointer, Process>(m,"ComputeNodalGradientProcess") + .def(py::init >& , Variable& >()) + .def(py::init&, Variable >& , Variable& >()) + ; + + m.attr("ComputeNodalGradientProcess2D") = m.attr("ComputeNodalGradientProcess"); + m.attr("ComputeNodalGradientProcess3D") = m.attr("ComputeNodalGradientProcess"); + m.attr("ComputeNodalGradientProcessComp2D") = m.attr("ComputeNodalGradientProcess"); + m.attr("ComputeNodalGradientProcessComp3D") = m.attr("ComputeNodalGradientProcess"); + /* Non-Historical */ - // DOUBLE - py::class_, NonHistorical>, ComputeNodalGradientProcess<2, Variable, NonHistorical>::Pointer, Process>(m,"ComputeNonHistoricalNodalGradientProcess2D") - .def(py::init&, Variable >& , Variable& >()) - ; - - py::class_, NonHistorical>, ComputeNodalGradientProcess<3, Variable, NonHistorical>::Pointer, Process>(m,"ComputeNonHistoricalNodalGradientProcess3D") - .def(py::init&, Variable >& , Variable& >()) - ; - - // COMPONENT - py::class_, ComputeNodalGradientProcess<2, component_type, NonHistorical>::Pointer, Process>(m,"ComputeNonHistoricalNodalGradientProcessComp2D") - .def(py::init >& , Variable& >()) + py::class_, ComputeNodalGradientProcess::Pointer, Process>(m,"ComputeNonHistoricalNodalGradientProcess") + .def(py::init >& , Variable& >()) + .def(py::init&, Variable >& , Variable& >()) ; - py::class_, ComputeNodalGradientProcess<3, component_type, NonHistorical>::Pointer, Process>(m,"ComputeNonHistoricalNodalGradientProcessComp3D") - .def(py::init >& , Variable& >()) - ; + m.attr("ComputeNonHistoricalNodalGradientProcess2D") = m.attr("ComputeNonHistoricalNodalGradientProcess"); + m.attr("ComputeNonHistoricalNodalGradientProcess3D") = m.attr("ComputeNonHistoricalNodalGradientProcess"); + m.attr("ComputeNonHistoricalNodalGradientProcessComp2D") = m.attr("ComputeNonHistoricalNodalGradientProcess"); + m.attr("ComputeNonHistoricalNodalGradientProcessComp3D") = m.attr("ComputeNonHistoricalNodalGradientProcess"); // Discontinuous distance computation methods py::class_, CalculateDiscontinuousDistanceToSkinProcess<2>::Pointer, Process>(m,"CalculateDiscontinuousDistanceToSkinProcess2D") diff --git a/kratos/tests/cpp_tests/processes/test_compute_nodal_gradient_process.cpp b/kratos/tests/cpp_tests/processes/test_compute_nodal_gradient_process.cpp index c997f530d98d..a49d3cbe7d90 100644 --- a/kratos/tests/cpp_tests/processes/test_compute_nodal_gradient_process.cpp +++ b/kratos/tests/cpp_tests/processes/test_compute_nodal_gradient_process.cpp @@ -17,7 +17,9 @@ // Project includes #include "containers/model.h" #include "geometries/triangle_2d_3.h" +#include "geometries/quadrilateral_2d_4.h" #include "geometries/tetrahedra_3d_4.h" +#include "geometries/hexahedra_3d_8.h" #include "testing/testing.h" #include "includes/kratos_flags.h" #include "includes/gid_io.h" @@ -45,7 +47,7 @@ namespace Kratos gid_io.WriteNodalResults(DISTANCE_GRADIENT, ThisModelPart.Nodes(), label, 0); } - /** + /** * Checks the correct work of the nodal gradient compute * Test triangle */ @@ -53,79 +55,84 @@ namespace Kratos KRATOS_TEST_CASE_IN_SUITE(NodalGradient1, KratosCoreFastSuite) { Model current_model; - + ModelPart& this_model_part = current_model.CreateModelPart("Main"); this_model_part.SetBufferSize(2); - + this_model_part.AddNodalSolutionStepVariable(DISTANCE); this_model_part.AddNodalSolutionStepVariable(DISTANCE_GRADIENT); - + Properties::Pointer p_elem_prop = this_model_part.pGetProperties(0); - + auto& process_info = this_model_part.GetProcessInfo(); - process_info[STEP] = 1; - process_info[NL_ITERATION_NUMBER] = 1; - - // First we create the nodes + process_info.SetValue(DOMAIN_SIZE, 2); + process_info.SetValue(STEP, 1); + process_info.SetValue(NL_ITERATION_NUMBER, 1); + + // First we create the nodes NodeType::Pointer p_node_1 = this_model_part.CreateNewNode(1, 0.0 , 0.0 , 0.0); NodeType::Pointer p_node_2 = this_model_part.CreateNewNode(2, 1.0 , 0.0 , 0.0); NodeType::Pointer p_node_3 = this_model_part.CreateNewNode(3, 1.0 , 1.0 , 0.0); NodeType::Pointer p_node_4 = this_model_part.CreateNewNode(4, 0.0 , 1.0 , 0.0); NodeType::Pointer p_node_5 = this_model_part.CreateNewNode(5, 2.0 , 0.0 , 0.0); NodeType::Pointer p_node_6 = this_model_part.CreateNewNode(6, 2.0 , 1.0 , 0.0); - + + // Initialize nodal area + for (auto& node : this_model_part.Nodes()) + node.SetValue(NODAL_AREA, 0.0); + // Now we create the "conditions" std::vector element_nodes_0 (3); element_nodes_0[0] = p_node_1; element_nodes_0[1] = p_node_2; element_nodes_0[2] = p_node_3; Triangle2D3 triangle_0( PointerVector{element_nodes_0} ); - + std::vector element_nodes_1 (3); element_nodes_1[0] = p_node_1; element_nodes_1[1] = p_node_3; element_nodes_1[2] = p_node_4; Triangle2D3 triangle_1( PointerVector{element_nodes_1} ); - + std::vector element_nodes_2 (3); element_nodes_2[0] = p_node_2; element_nodes_2[1] = p_node_5; element_nodes_2[2] = p_node_3; Triangle2D3 triangle_2( PointerVector{element_nodes_2} ); - + std::vector element_nodes_3 (3); element_nodes_3[0] = p_node_5; element_nodes_3[1] = p_node_6; element_nodes_3[2] = p_node_3; Triangle2D3 triangle_3( PointerVector{element_nodes_3} ); - + Element::Pointer p_elem_0 = this_model_part.CreateNewElement("Element2D3N", 1, triangle_0, p_elem_prop); Element::Pointer p_elem_1 = this_model_part.CreateNewElement("Element2D3N", 2, triangle_1, p_elem_prop); Element::Pointer p_elem_2 = this_model_part.CreateNewElement("Element2D3N", 3, triangle_2, p_elem_prop); Element::Pointer p_elem_3 = this_model_part.CreateNewElement("Element2D3N", 4, triangle_3, p_elem_prop); - + // Set DISTANCE for (std::size_t i_node = 0; i_node < this_model_part.Nodes().size(); ++i_node) { 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_AREA, 0.0); } - - typedef ComputeNodalGradientProcess<2, Variable, Historical> GradientType; + + typedef ComputeNodalGradientProcess GradientType; GradientType process = GradientType(this_model_part, DISTANCE, DISTANCE_GRADIENT); process.Execute(); - -// // DEBUG + +// // DEBUG // GiDIODebugGradient(this_model_part); - + const double tolerance = 1.0e-8; KRATOS_CHECK_LESS_EQUAL(std::abs(p_node_1->FastGetSolutionStepValue(DISTANCE_GRADIENT_X)) - 1.0, tolerance); KRATOS_CHECK_LESS_EQUAL(std::abs(p_node_2->FastGetSolutionStepValue(DISTANCE_GRADIENT_X)) - 1.0, tolerance); KRATOS_CHECK_LESS_EQUAL(std::abs(p_node_5->FastGetSolutionStepValue(DISTANCE_GRADIENT_X)) - 1.0, tolerance); KRATOS_CHECK_LESS_EQUAL(std::abs(p_node_6->FastGetSolutionStepValue(DISTANCE_GRADIENT_X)) - 1.0, tolerance); } - - /** + + /** * Checks the correct work of the nodal gradient compute * Test tetrahedra */ @@ -136,31 +143,36 @@ namespace Kratos ModelPart& this_model_part = current_model.CreateModelPart("Main"); this_model_part.SetBufferSize(2); - + this_model_part.AddNodalSolutionStepVariable(DISTANCE); this_model_part.AddNodalSolutionStepVariable(DISTANCE_GRADIENT); - + Properties::Pointer p_elem_prop = this_model_part.pGetProperties(0); - + auto& process_info = this_model_part.GetProcessInfo(); - process_info[STEP] = 1; - process_info[NL_ITERATION_NUMBER] = 1; - - // First we create the nodes + process_info.SetValue(DOMAIN_SIZE, 3); + process_info.SetValue(STEP, 1); + process_info.SetValue(NL_ITERATION_NUMBER, 1); + + // First we create the nodes NodeType::Pointer p_node_1 = this_model_part.CreateNewNode(1 , 0.0 , 1.0 , 1.0); NodeType::Pointer p_node_2 = this_model_part.CreateNewNode(2 , 0.0 , 1.0 , 0.0); NodeType::Pointer p_node_3 = this_model_part.CreateNewNode(3 , 0.0 , 0.0 , 1.0); NodeType::Pointer p_node_4 = this_model_part.CreateNewNode(4 , 1.0 , 1.0 , 1.0); NodeType::Pointer p_node_5 = this_model_part.CreateNewNode(5 , 0.0 , 0.0 , 0.0); NodeType::Pointer p_node_6 = this_model_part.CreateNewNode(6 , 1.0 , 1.0 , 0.0); - + NodeType::Pointer p_node_7 = this_model_part.CreateNewNode(7 , 1.0 , 0.0 , 1.0); NodeType::Pointer p_node_8 = this_model_part.CreateNewNode(8 , 1.0 , 0.0 , 0.0); NodeType::Pointer p_node_9 = this_model_part.CreateNewNode(9 , 2.0 , 1.0 , 1.0); NodeType::Pointer p_node_10 = this_model_part.CreateNewNode(10 , 2.0 , 1.0 , 0.0); NodeType::Pointer p_node_11 = this_model_part.CreateNewNode(11 , 2.0 , 0.0 , 1.0); NodeType::Pointer p_node_12 = this_model_part.CreateNewNode(12 , 2.0 , 0.0 , 0.0); - + + // Initialize nodal area + for (auto& node : this_model_part.Nodes()) + node.SetValue(NODAL_AREA, 0.0); + // Now we create the "conditions" std::vector element_nodes_0 (4); element_nodes_0[0] = p_node_12; @@ -168,84 +180,84 @@ namespace Kratos element_nodes_0[2] = p_node_8; element_nodes_0[3] = p_node_9; Tetrahedra3D4 tetrahedra_0( PointerVector{element_nodes_0} ); - + std::vector element_nodes_1 (4); element_nodes_1[0] = p_node_4; element_nodes_1[1] = p_node_6; element_nodes_1[2] = p_node_9; element_nodes_1[3] = p_node_7; Tetrahedra3D4 tetrahedra_1( PointerVector{element_nodes_1} ); - + std::vector element_nodes_2 (4); element_nodes_2[0] = p_node_11; element_nodes_2[1] = p_node_7; element_nodes_2[2] = p_node_9; element_nodes_2[3] = p_node_8; Tetrahedra3D4 tetrahedra_2( PointerVector{element_nodes_2} ); - + std::vector element_nodes_3 (4); element_nodes_3[0] = p_node_5; element_nodes_3[1] = p_node_3; element_nodes_3[2] = p_node_8; element_nodes_3[3] = p_node_6; Tetrahedra3D4 tetrahedra_3( PointerVector{element_nodes_3} ); - + std::vector element_nodes_4 (4); element_nodes_4[0] = p_node_4; element_nodes_4[1] = p_node_6; element_nodes_4[2] = p_node_7; element_nodes_4[3] = p_node_3; Tetrahedra3D4 tetrahedra_4( PointerVector{element_nodes_4} ); - + std::vector element_nodes_5 (4); element_nodes_5[0] = p_node_2; element_nodes_5[1] = p_node_3; element_nodes_5[2] = p_node_5; element_nodes_5[3] = p_node_6; Tetrahedra3D4 tetrahedra_5( PointerVector{element_nodes_5} ); - + std::vector element_nodes_6 (4); element_nodes_6[0] = p_node_10; element_nodes_6[1] = p_node_9; element_nodes_6[2] = p_node_6; element_nodes_6[3] = p_node_8; Tetrahedra3D4 tetrahedra_6( PointerVector{element_nodes_6} ); - + std::vector element_nodes_7 (4); element_nodes_7[0] = p_node_7; element_nodes_7[1] = p_node_8; element_nodes_7[2] = p_node_3; element_nodes_7[3] = p_node_6; Tetrahedra3D4 tetrahedra_7( PointerVector{element_nodes_7} ); - + std::vector element_nodes_8 (4); element_nodes_8[0] = p_node_7; element_nodes_8[1] = p_node_8; element_nodes_8[2] = p_node_6; element_nodes_8[3] = p_node_9; Tetrahedra3D4 tetrahedra_8( PointerVector{element_nodes_8} ); - + std::vector element_nodes_9 (4); element_nodes_9[0] = p_node_4; element_nodes_9[1] = p_node_1; element_nodes_9[2] = p_node_6; element_nodes_9[3] = p_node_3; Tetrahedra3D4 tetrahedra_9( PointerVector{element_nodes_9} ); - + std::vector element_nodes_10 (4); element_nodes_10[0] = p_node_9; element_nodes_10[1] = p_node_12; element_nodes_10[2] = p_node_11; element_nodes_10[3] = p_node_8; Tetrahedra3D4 tetrahedra_10( PointerVector{element_nodes_10} ); - + std::vector element_nodes_11 (4); element_nodes_11[0] = p_node_3; element_nodes_11[1] = p_node_2; element_nodes_11[2] = p_node_1; element_nodes_11[3] = p_node_6; Tetrahedra3D4 tetrahedra_11( PointerVector{element_nodes_11} ); - + Element::Pointer p_elem_0 = this_model_part.CreateNewElement("Element3D4N", 1, tetrahedra_0, p_elem_prop); Element::Pointer p_elem_1 = this_model_part.CreateNewElement("Element3D4N", 2, tetrahedra_1, p_elem_prop); Element::Pointer p_elem_2 = this_model_part.CreateNewElement("Element3D4N", 3, tetrahedra_2, p_elem_prop); @@ -258,22 +270,185 @@ namespace Kratos Element::Pointer p_elem_9 = this_model_part.CreateNewElement("Element3D4N", 10, tetrahedra_9, p_elem_prop); Element::Pointer p_elem_10 = this_model_part.CreateNewElement("Element3D4N", 11, tetrahedra_10, p_elem_prop); Element::Pointer p_elem_11 = this_model_part.CreateNewElement("Element3D4N", 12, tetrahedra_11, p_elem_prop); - + // Set DISTANCE for (std::size_t i_node = 0; i_node < this_model_part.Nodes().size(); ++i_node) { 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_AREA, 0.0); } - + // Compute gradient - typedef ComputeNodalGradientProcess<3, Variable, Historical> GradientType; + typedef ComputeNodalGradientProcess GradientType; + GradientType process = GradientType(this_model_part, DISTANCE, DISTANCE_GRADIENT); + process.Execute(); + +// // DEBUG +// GiDIODebugGradient(this_model_part); + + const double tolerance = 1.0e-8; + KRATOS_CHECK_LESS_EQUAL(std::abs(p_node_1->FastGetSolutionStepValue(DISTANCE_GRADIENT_X)) - 1.0, tolerance); + KRATOS_CHECK_LESS_EQUAL(std::abs(p_node_2->FastGetSolutionStepValue(DISTANCE_GRADIENT_X)) - 1.0, tolerance); + KRATOS_CHECK_LESS_EQUAL(std::abs(p_node_3->FastGetSolutionStepValue(DISTANCE_GRADIENT_X)) - 1.0, tolerance); + KRATOS_CHECK_LESS_EQUAL(std::abs(p_node_5->FastGetSolutionStepValue(DISTANCE_GRADIENT_X)) - 1.0, tolerance); + KRATOS_CHECK_LESS_EQUAL(std::abs(p_node_9->FastGetSolutionStepValue(DISTANCE_GRADIENT_X)) - 1.0, tolerance); + KRATOS_CHECK_LESS_EQUAL(std::abs(p_node_10->FastGetSolutionStepValue(DISTANCE_GRADIENT_X)) - 1.0, tolerance); + KRATOS_CHECK_LESS_EQUAL(std::abs(p_node_11->FastGetSolutionStepValue(DISTANCE_GRADIENT_X)) - 1.0, tolerance); + KRATOS_CHECK_LESS_EQUAL(std::abs(p_node_12->FastGetSolutionStepValue(DISTANCE_GRADIENT_X)) - 1.0, tolerance); + } + + /** + * Checks the correct work of the nodal gradient compute + * Test quadrilateral + */ + KRATOS_TEST_CASE_IN_SUITE(NodalGradient3, KratosCoreFastSuite) + { + Model current_model; + + ModelPart& this_model_part = current_model.CreateModelPart("Main"); + this_model_part.SetBufferSize(2); + + this_model_part.AddNodalSolutionStepVariable(DISTANCE); + this_model_part.AddNodalSolutionStepVariable(DISTANCE_GRADIENT); + + Properties::Pointer p_elem_prop = this_model_part.pGetProperties(0); + + auto& process_info = this_model_part.GetProcessInfo(); + process_info.SetValue(DOMAIN_SIZE, 2); + process_info.SetValue(STEP, 1); + process_info.SetValue(NL_ITERATION_NUMBER, 1); + + // First we create the nodes + NodeType::Pointer p_node_1 = this_model_part.CreateNewNode(1, 0.0 , 0.0 , 0.0); + NodeType::Pointer p_node_2 = this_model_part.CreateNewNode(2, 1.0 , 0.0 , 0.0); + NodeType::Pointer p_node_3 = this_model_part.CreateNewNode(3, 1.0 , 1.0 , 0.0); + NodeType::Pointer p_node_4 = this_model_part.CreateNewNode(4, 0.0 , 1.0 , 0.0); + NodeType::Pointer p_node_5 = this_model_part.CreateNewNode(5, 2.0 , 0.0 , 0.0); + NodeType::Pointer p_node_6 = this_model_part.CreateNewNode(6, 2.0 , 1.0 , 0.0); + + // Initialize nodal area + for (auto& node : this_model_part.Nodes()) + node.SetValue(NODAL_AREA, 0.0); + + // Now we create the "conditions" + std::vector element_nodes_0 (4); + element_nodes_0[0] = p_node_1; + element_nodes_0[1] = p_node_2; + element_nodes_0[2] = p_node_3; + element_nodes_0[3] = p_node_4; + Quadrilateral2D4 quadrilateral_0( PointerVector{element_nodes_0} ); + + std::vector element_nodes_1 (4); + element_nodes_1[0] = p_node_2; + element_nodes_1[1] = p_node_5; + element_nodes_1[2] = p_node_6; + element_nodes_1[3] = p_node_3; + Quadrilateral2D4 quadrilateral_1( PointerVector{element_nodes_1} ); + + Element::Pointer p_elem_0 = this_model_part.CreateNewElement("Element2D4N", 1, quadrilateral_0, p_elem_prop); + Element::Pointer p_elem_1 = this_model_part.CreateNewElement("Element2D4N", 2, quadrilateral_1, p_elem_prop); + + // Set DISTANCE + for (std::size_t i_node = 0; i_node < this_model_part.Nodes().size(); ++i_node) { + 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_AREA, 0.0); + } + + typedef ComputeNodalGradientProcess GradientType; GradientType process = GradientType(this_model_part, DISTANCE, DISTANCE_GRADIENT); process.Execute(); // // DEBUG // GiDIODebugGradient(this_model_part); + const double tolerance = 1.0e-8; + KRATOS_CHECK_LESS_EQUAL(std::abs(p_node_1->FastGetSolutionStepValue(DISTANCE_GRADIENT_X)) - 1.0, tolerance); + KRATOS_CHECK_LESS_EQUAL(std::abs(p_node_2->FastGetSolutionStepValue(DISTANCE_GRADIENT_X)) - 1.0, tolerance); + KRATOS_CHECK_LESS_EQUAL(std::abs(p_node_5->FastGetSolutionStepValue(DISTANCE_GRADIENT_X)) - 1.0, tolerance); + KRATOS_CHECK_LESS_EQUAL(std::abs(p_node_6->FastGetSolutionStepValue(DISTANCE_GRADIENT_X)) - 1.0, tolerance); + } + + /** + * Checks the correct work of the nodal gradient compute + * Test hexahedra + */ + KRATOS_TEST_CASE_IN_SUITE(NodalGradient4, KratosCoreFastSuite) + { + Model current_model; + + ModelPart& this_model_part = current_model.CreateModelPart("Main"); + this_model_part.SetBufferSize(2); + + this_model_part.AddNodalSolutionStepVariable(DISTANCE); + this_model_part.AddNodalSolutionStepVariable(DISTANCE_GRADIENT); + + Properties::Pointer p_elem_prop = this_model_part.pGetProperties(0); + + auto& process_info = this_model_part.GetProcessInfo(); + process_info.SetValue(DOMAIN_SIZE, 3); + process_info.SetValue(STEP, 1); + process_info.SetValue(NL_ITERATION_NUMBER, 1); + + // First we create the nodes + NodeType::Pointer p_node_1 = this_model_part.CreateNewNode(1 , 0.0 , 1.0 , 1.0); + NodeType::Pointer p_node_2 = this_model_part.CreateNewNode(2 , 0.0 , 1.0 , 0.0); + NodeType::Pointer p_node_3 = this_model_part.CreateNewNode(3 , 0.0 , 0.0 , 1.0); + NodeType::Pointer p_node_4 = this_model_part.CreateNewNode(4 , 1.0 , 1.0 , 1.0); + NodeType::Pointer p_node_5 = this_model_part.CreateNewNode(5 , 0.0 , 0.0 , 0.0); + NodeType::Pointer p_node_6 = this_model_part.CreateNewNode(6 , 1.0 , 1.0 , 0.0); + NodeType::Pointer p_node_7 = this_model_part.CreateNewNode(7 , 1.0 , 0.0 , 1.0); + NodeType::Pointer p_node_8 = this_model_part.CreateNewNode(8 , 1.0 , 0.0 , 0.0); + NodeType::Pointer p_node_9 = this_model_part.CreateNewNode(9 , 2.0 , 1.0 , 1.0); + NodeType::Pointer p_node_10 = this_model_part.CreateNewNode(10 , 2.0 , 1.0 , 0.0); + NodeType::Pointer p_node_11 = this_model_part.CreateNewNode(11 , 2.0 , 0.0 , 1.0); + NodeType::Pointer p_node_12 = this_model_part.CreateNewNode(12 , 2.0 , 0.0 , 0.0); + + // Initialize nodal area + for (auto& node : this_model_part.Nodes()) + node.SetValue(NODAL_AREA, 0.0); + + // Now we create the "conditions" + std::vector element_nodes_0 (8); + element_nodes_0[0] = p_node_5; + element_nodes_0[1] = p_node_8; + element_nodes_0[2] = p_node_6; + element_nodes_0[3] = p_node_2; + element_nodes_0[4] = p_node_3; + element_nodes_0[5] = p_node_7; + element_nodes_0[6] = p_node_4; + element_nodes_0[7] = p_node_1; + Hexahedra3D8 hexahedra_0( PointerVector{element_nodes_0} ); + + std::vector element_nodes_1 (8); + element_nodes_1[0] = p_node_8; + element_nodes_1[1] = p_node_12; + element_nodes_1[2] = p_node_10; + element_nodes_1[3] = p_node_6; + element_nodes_1[4] = p_node_7; + element_nodes_1[5] = p_node_11; + element_nodes_1[6] = p_node_9; + element_nodes_1[7] = p_node_4; + Hexahedra3D8 hexahedra_1( PointerVector{element_nodes_1} ); + + Element::Pointer p_elem_0 = this_model_part.CreateNewElement("Element3D8N", 1, hexahedra_0, p_elem_prop); + Element::Pointer p_elem_1 = this_model_part.CreateNewElement("Element3D8N", 2, hexahedra_1, p_elem_prop); + + // Set DISTANCE + for (std::size_t i_node = 0; i_node < this_model_part.Nodes().size(); ++i_node) { + 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_AREA, 0.0); + } + + // Compute gradient + typedef ComputeNodalGradientProcess GradientType; + GradientType process = GradientType(this_model_part, DISTANCE, DISTANCE_GRADIENT); + process.Execute(); + +// // DEBUG +// GiDIODebugGradient(this_model_part); + const double tolerance = 1.0e-8; KRATOS_CHECK_LESS_EQUAL(std::abs(p_node_1->FastGetSolutionStepValue(DISTANCE_GRADIENT_X)) - 1.0, tolerance); KRATOS_CHECK_LESS_EQUAL(std::abs(p_node_2->FastGetSolutionStepValue(DISTANCE_GRADIENT_X)) - 1.0, tolerance);