diff --git a/applications/StructuralMechanicsApplication/custom_response_functions/adjoint_conditions/adjoint_semi_analytic_base_condition.cpp b/applications/StructuralMechanicsApplication/custom_response_functions/adjoint_conditions/adjoint_semi_analytic_base_condition.cpp index 6391f02032f9..a764dac0d8c2 100644 --- a/applications/StructuralMechanicsApplication/custom_response_functions/adjoint_conditions/adjoint_semi_analytic_base_condition.cpp +++ b/applications/StructuralMechanicsApplication/custom_response_functions/adjoint_conditions/adjoint_semi_analytic_base_condition.cpp @@ -28,6 +28,40 @@ namespace Kratos { + namespace AdjointSemiAnalyticBaseConditionHelperUtils + { + template + void CalculateOnIntegrationPoints( + Condition& rPrimalCondition, + const Condition& rAdjointCondition, + const Variable& rVariable, + std::vector& rValues, + const ProcessInfo& rCurrentProcessInfo) + { + KRATOS_TRY + + if (rAdjointCondition.Has(rVariable)) { + // Get result value for output + const auto& output_value = rAdjointCondition.GetValue(rVariable); + + // Resize Output + const auto gauss_points_number = rAdjointCondition.GetGeometry().IntegrationPointsNumber(rAdjointCondition.GetIntegrationMethod()); + if (rValues.size() != gauss_points_number) { + rValues.resize(gauss_points_number); + } + + // Write scalar result value on all Gauss-Points + for (IndexType i = 0; i < gauss_points_number; ++i) { + rValues[i] = output_value; + } + } + else { + rPrimalCondition.CalculateOnIntegrationPoints(rVariable, rValues, rCurrentProcessInfo); + } + + KRATOS_CATCH(""); + } + } // namespace AdjointSemiAnalyticBaseConditionHelperUtils template void AdjointSemiAnalyticBaseCondition::EquationIdVector(EquationIdVectorType& rResult, const ProcessInfo& rCurrentProcessInfo ) const @@ -114,62 +148,74 @@ namespace Kratos template void AdjointSemiAnalyticBaseCondition::CalculateOnIntegrationPoints( - const Variable >& rVariable, std::vector< array_1d >& rOutput, const ProcessInfo& rCurrentProcessInfo) + const Variable& rVariable, + std::vector& rOutput, + const ProcessInfo& rCurrentProcessInfo) { - KRATOS_TRY; - - if (this->Has(rVariable)) { - // Get result value for output - const auto& output_value = this->GetValue(rVariable); - - // Resize Output - const SizeType gauss_points_number = this->GetGeometry() - .IntegrationPointsNumber(this->GetIntegrationMethod()); - if (rOutput.size() != gauss_points_number) { - rOutput.resize(gauss_points_number); - } - - // Write result value on all Gauss-Points - for(IndexType i = 0; i < gauss_points_number; ++i) { - rOutput[i] = output_value; - } - - } - else { - KRATOS_ERROR << "Unsupported output variable." << std::endl; - } + AdjointSemiAnalyticBaseConditionHelperUtils::CalculateOnIntegrationPoints(*mpPrimalCondition, *this, rVariable, rOutput, rCurrentProcessInfo); + } - KRATOS_CATCH("") + template + void AdjointSemiAnalyticBaseCondition::CalculateOnIntegrationPoints( + const Variable& rVariable, + std::vector& rOutput, + const ProcessInfo& rCurrentProcessInfo) + { + AdjointSemiAnalyticBaseConditionHelperUtils::CalculateOnIntegrationPoints(*mpPrimalCondition, *this, rVariable, rOutput, rCurrentProcessInfo); } template void AdjointSemiAnalyticBaseCondition::CalculateOnIntegrationPoints( - const Variable& rVariable, std::vector& rOutput, const ProcessInfo& rCurrentProcessInfo) + const Variable>& rVariable, + std::vector>& rOutput, + const ProcessInfo& rCurrentProcessInfo) { - KRATOS_TRY; + AdjointSemiAnalyticBaseConditionHelperUtils::CalculateOnIntegrationPoints(*mpPrimalCondition, *this, rVariable, rOutput, rCurrentProcessInfo); + } - if (this->Has(rVariable)) { - // Get result value for output - const auto& output_value = this->GetValue(rVariable); + template + void AdjointSemiAnalyticBaseCondition::CalculateOnIntegrationPoints( + const Variable>& rVariable, + std::vector>& rOutput, + const ProcessInfo& rCurrentProcessInfo) + { + AdjointSemiAnalyticBaseConditionHelperUtils::CalculateOnIntegrationPoints(*mpPrimalCondition, *this, rVariable, rOutput, rCurrentProcessInfo); + } - // Resize Output - const SizeType gauss_points_number = this->GetGeometry() - .IntegrationPointsNumber(this->GetIntegrationMethod()); - if (rOutput.size() != gauss_points_number) { - rOutput.resize(gauss_points_number); - } + template + void AdjointSemiAnalyticBaseCondition::CalculateOnIntegrationPoints( + const Variable>& rVariable, + std::vector>& rOutput, + const ProcessInfo& rCurrentProcessInfo) + { + AdjointSemiAnalyticBaseConditionHelperUtils::CalculateOnIntegrationPoints(*mpPrimalCondition, *this, rVariable, rOutput, rCurrentProcessInfo); + } - // Write result value on all Gauss-Points - for(IndexType i = 0; i < gauss_points_number; ++i) { - rOutput[i] = output_value; - } + template + void AdjointSemiAnalyticBaseCondition::CalculateOnIntegrationPoints( + const Variable>& rVariable, + std::vector>& rOutput, + const ProcessInfo& rCurrentProcessInfo) + { + AdjointSemiAnalyticBaseConditionHelperUtils::CalculateOnIntegrationPoints(*mpPrimalCondition, *this, rVariable, rOutput, rCurrentProcessInfo); + } - } - else { - KRATOS_ERROR << "Unsupported output variable." << std::endl; - } + template + void AdjointSemiAnalyticBaseCondition::CalculateOnIntegrationPoints( + const Variable& rVariable, + std::vector& rOutput, + const ProcessInfo& rCurrentProcessInfo) + { + AdjointSemiAnalyticBaseConditionHelperUtils::CalculateOnIntegrationPoints(*mpPrimalCondition, *this, rVariable, rOutput, rCurrentProcessInfo); + } - KRATOS_CATCH("") + template + void AdjointSemiAnalyticBaseCondition::CalculateOnIntegrationPoints( + const Variable& rVariable, + std::vector& rOutput, + const ProcessInfo& rCurrentProcessInfo) + { + AdjointSemiAnalyticBaseConditionHelperUtils::CalculateOnIntegrationPoints(*mpPrimalCondition, *this, rVariable, rOutput, rCurrentProcessInfo); } template diff --git a/applications/StructuralMechanicsApplication/custom_response_functions/adjoint_conditions/adjoint_semi_analytic_base_condition.h b/applications/StructuralMechanicsApplication/custom_response_functions/adjoint_conditions/adjoint_semi_analytic_base_condition.h index a44da912cf97..8d38e0028ded 100644 --- a/applications/StructuralMechanicsApplication/custom_response_functions/adjoint_conditions/adjoint_semi_analytic_base_condition.h +++ b/applications/StructuralMechanicsApplication/custom_response_functions/adjoint_conditions/adjoint_semi_analytic_base_condition.h @@ -258,28 +258,46 @@ class AdjointSemiAnalyticBaseCondition KRATOS_ERROR << "Calculate of the adjoint base condition is called!" << std::endl; } - - void CalculateOnIntegrationPoints(const Variable& rVariable, - std::vector& rOutput, - const ProcessInfo& rCurrentProcessInfo) override; - - void CalculateOnIntegrationPoints(const Variable >& rVariable, - std::vector< array_1d >& Output, - const ProcessInfo& rCurrentProcessInfo) override; - - void CalculateOnIntegrationPoints(const Variable& rVariable, - std::vector< Vector >& Output, - const ProcessInfo& rCurrentProcessInfo) override - { - KRATOS_ERROR << "CalculateOnIntegrationPoints of the adjoint base condition is called!" << std::endl; - } - - void CalculateOnIntegrationPoints(const Variable& rVariable, - std::vector< Matrix >& Output, - const ProcessInfo& rCurrentProcessInfo) override - { - KRATOS_ERROR << "CalculateOnIntegrationPoints of the adjoint base condition is called!" << std::endl; - } + // Results calculation on integration points + void CalculateOnIntegrationPoints( + const Variable& rVariable, + std::vector& rOutput, + const ProcessInfo& rCurrentProcessInfo) override; + + void CalculateOnIntegrationPoints( + const Variable& rVariable, + std::vector& rOutput, + const ProcessInfo& rCurrentProcessInfo) override; + + void CalculateOnIntegrationPoints( + const Variable>& rVariable, + std::vector>& rOutput, + const ProcessInfo& rCurrentProcessInfo) override; + + void CalculateOnIntegrationPoints( + const Variable>& rVariable, + std::vector>& rOutput, + const ProcessInfo& rCurrentProcessInfo) override; + + void CalculateOnIntegrationPoints( + const Variable>& rVariable, + std::vector>& rOutput, + const ProcessInfo& rCurrentProcessInfo) override; + + void CalculateOnIntegrationPoints( + const Variable>& rVariable, + std::vector>& rOutput, + const ProcessInfo& rCurrentProcessInfo) override; + + void CalculateOnIntegrationPoints( + const Variable& rVariable, + std::vector& rOutput, + const ProcessInfo& rCurrentProcessInfo) override; + + void CalculateOnIntegrationPoints( + const Variable& rVariable, + std::vector& rOutput, + const ProcessInfo& rCurrentProcessInfo) override; int Check( const ProcessInfo& rCurrentProcessInfo ) const override; diff --git a/applications/StructuralMechanicsApplication/custom_response_functions/adjoint_elements/adjoint_finite_difference_base_element.cpp b/applications/StructuralMechanicsApplication/custom_response_functions/adjoint_elements/adjoint_finite_difference_base_element.cpp index e1b8bda6b227..89cb7d2ecce8 100644 --- a/applications/StructuralMechanicsApplication/custom_response_functions/adjoint_elements/adjoint_finite_difference_base_element.cpp +++ b/applications/StructuralMechanicsApplication/custom_response_functions/adjoint_elements/adjoint_finite_difference_base_element.cpp @@ -29,6 +29,42 @@ namespace Kratos { +namespace AdjointFiniteDifferenceBaseElementHelperUtils +{ + +template +void CalculateOnIntegrationPoints( + Element& rPrimalElement, + const Element& rAdjointElement, + const Variable& rVariable, + std::vector& rValues, + const ProcessInfo& rCurrentProcessInfo) +{ + KRATOS_TRY + + if (rAdjointElement.Has(rVariable)) { + // Get result value for output + const auto& output_value = rAdjointElement.GetValue(rVariable); + + // Resize Output + const SizeType gauss_points_number = rAdjointElement.GetGeometry().IntegrationPointsNumber(rAdjointElement.GetIntegrationMethod()); + if (rValues.size() != gauss_points_number) { + rValues.resize(gauss_points_number); + } + + // Write scalar result value on all Gauss-Points + for (IndexType i = 0; i < gauss_points_number; ++i) { + rValues[i] = output_value; + } + } + else { + rPrimalElement.CalculateOnIntegrationPoints(rVariable, rValues, rCurrentProcessInfo); + } + + KRATOS_CATCH(""); +} +} // namespace AdjointFiniteDifferenceBaseElementHelperUtils + template void AdjointFiniteDifferencingBaseElement::EquationIdVector(EquationIdVectorType& rResult, const ProcessInfo& rCurrentProcessInfo) const @@ -194,60 +230,75 @@ void AdjointFiniteDifferencingBaseElement::Calculate(const Varia } template -void AdjointFiniteDifferencingBaseElement::CalculateOnIntegrationPoints(const Variable& rVariable, - std::vector& rValues, - const ProcessInfo& rCurrentProcessInfo) +void AdjointFiniteDifferencingBaseElement::CalculateOnIntegrationPoints( + const Variable& rVariable, + std::vector& rOutput, + const ProcessInfo& rCurrentProcessInfo) { - KRATOS_TRY; - - if(this->Has(rVariable)) - { - // Get result value for output - const double& output_value = this->GetValue(rVariable); - - // Resize Output - const SizeType gauss_points_number = this->GetGeometry() - .IntegrationPointsNumber(this->GetIntegrationMethod()); - if (rValues.size() != gauss_points_number) - rValues.resize(gauss_points_number); - - // Write scalar result value on all Gauss-Points - for(IndexType i = 0; i < gauss_points_number; ++i) - rValues[i] = output_value; - } - else - KRATOS_ERROR << "Unsupported output variable." << std::endl; + AdjointFiniteDifferenceBaseElementHelperUtils::CalculateOnIntegrationPoints(*mpPrimalElement, *this, rVariable, rOutput, rCurrentProcessInfo); +} - KRATOS_CATCH("") +template +void AdjointFiniteDifferencingBaseElement::CalculateOnIntegrationPoints( + const Variable& rVariable, + std::vector& rOutput, + const ProcessInfo& rCurrentProcessInfo) +{ + AdjointFiniteDifferenceBaseElementHelperUtils::CalculateOnIntegrationPoints(*mpPrimalElement, *this, rVariable, rOutput, rCurrentProcessInfo); } template void AdjointFiniteDifferencingBaseElement::CalculateOnIntegrationPoints( - const Variable >& rVariable, std::vector< array_1d >& rOutput, const ProcessInfo& rCurrentProcessInfo) + const Variable>& rVariable, + std::vector>& rOutput, + const ProcessInfo& rCurrentProcessInfo) { - KRATOS_TRY; + AdjointFiniteDifferenceBaseElementHelperUtils::CalculateOnIntegrationPoints(*mpPrimalElement, *this, rVariable, rOutput, rCurrentProcessInfo); +} - if(this->Has(rVariable)) { - // Get result value for output - const auto& output_value = this->GetValue(rVariable); +template +void AdjointFiniteDifferencingBaseElement::CalculateOnIntegrationPoints( + const Variable>& rVariable, + std::vector>& rOutput, + const ProcessInfo& rCurrentProcessInfo) +{ + AdjointFiniteDifferenceBaseElementHelperUtils::CalculateOnIntegrationPoints(*mpPrimalElement, *this, rVariable, rOutput, rCurrentProcessInfo); +} - // Resize Output - const SizeType gauss_points_number = this->GetGeometry() - .IntegrationPointsNumber(this->GetIntegrationMethod()); - if (rOutput.size() != gauss_points_number) { - rOutput.resize(gauss_points_number); - } +template +void AdjointFiniteDifferencingBaseElement::CalculateOnIntegrationPoints( + const Variable>& rVariable, + std::vector>& rOutput, + const ProcessInfo& rCurrentProcessInfo) +{ + AdjointFiniteDifferenceBaseElementHelperUtils::CalculateOnIntegrationPoints(*mpPrimalElement, *this, rVariable, rOutput, rCurrentProcessInfo); +} - // Write scalar result value on all Gauss-Points - for(IndexType i = 0; i < gauss_points_number; ++i) { - rOutput[i] = output_value; - } +template +void AdjointFiniteDifferencingBaseElement::CalculateOnIntegrationPoints( + const Variable>& rVariable, + std::vector>& rOutput, + const ProcessInfo& rCurrentProcessInfo) +{ + AdjointFiniteDifferenceBaseElementHelperUtils::CalculateOnIntegrationPoints(*mpPrimalElement, *this, rVariable, rOutput, rCurrentProcessInfo); +} - } else { - KRATOS_ERROR << "Unsupported output variable." << std::endl; - } +template +void AdjointFiniteDifferencingBaseElement::CalculateOnIntegrationPoints( + const Variable& rVariable, + std::vector& rOutput, + const ProcessInfo& rCurrentProcessInfo) +{ + AdjointFiniteDifferenceBaseElementHelperUtils::CalculateOnIntegrationPoints(*mpPrimalElement, *this, rVariable, rOutput, rCurrentProcessInfo); +} - KRATOS_CATCH("") +template +void AdjointFiniteDifferencingBaseElement::CalculateOnIntegrationPoints( + const Variable& rVariable, + std::vector& rOutput, + const ProcessInfo& rCurrentProcessInfo) +{ + AdjointFiniteDifferenceBaseElementHelperUtils::CalculateOnIntegrationPoints(*mpPrimalElement, *this, rVariable, rOutput, rCurrentProcessInfo); } template diff --git a/applications/StructuralMechanicsApplication/custom_response_functions/adjoint_elements/adjoint_finite_difference_base_element.h b/applications/StructuralMechanicsApplication/custom_response_functions/adjoint_elements/adjoint_finite_difference_base_element.h index 83b7b9c64860..2a1da05de185 100644 --- a/applications/StructuralMechanicsApplication/custom_response_functions/adjoint_elements/adjoint_finite_difference_base_element.h +++ b/applications/StructuralMechanicsApplication/custom_response_functions/adjoint_elements/adjoint_finite_difference_base_element.h @@ -304,42 +304,45 @@ class AdjointFiniteDifferencingBaseElement : public Element const ProcessInfo& rCurrentProcessInfo) override; // Results calculation on integration points - void CalculateOnIntegrationPoints(const Variable& rVariable, - std::vector& rOutput, - const ProcessInfo& rCurrentProcessInfo) override - { - KRATOS_ERROR << "CalculateOnIntegrationPoints of the adjoint base element is called!" << std::endl; - } - - void CalculateOnIntegrationPoints(const Variable& rVariable, - std::vector& rOutput, - const ProcessInfo& rCurrentProcessInfo) override; - - void CalculateOnIntegrationPoints(const Variable >& rVariable, - std::vector< array_1d >& rOutput, - const ProcessInfo& rCurrentProcessInfo) override; - - void CalculateOnIntegrationPoints(const Variable >& rVariable, - std::vector< array_1d >& rOutput, - const ProcessInfo& rCurrentProcessInfo) override - { - KRATOS_ERROR << "CalculateOnIntegrationPoints of the adjoint base element is called!" << std::endl; - } - - void CalculateOnIntegrationPoints(const Variable& rVariable, - std::vector< Vector >& rOutput, - const ProcessInfo& rCurrentProcessInfo) override - { - KRATOS_ERROR << "CalculateOnIntegrationPoints of the adjoint base element is called!" << std::endl; - } - - void CalculateOnIntegrationPoints(const Variable& rVariable, - std::vector< Matrix >& rOutput, - const ProcessInfo& rCurrentProcessInfo) override - { - KRATOS_ERROR << "CalculateOnIntegrationPoints of the adjoint base element is called!" << std::endl; - } - + void CalculateOnIntegrationPoints( + const Variable& rVariable, + std::vector& rOutput, + const ProcessInfo& rCurrentProcessInfo) override; + + void CalculateOnIntegrationPoints( + const Variable& rVariable, + std::vector& rOutput, + const ProcessInfo& rCurrentProcessInfo) override; + + void CalculateOnIntegrationPoints( + const Variable>& rVariable, + std::vector>& rOutput, + const ProcessInfo& rCurrentProcessInfo) override; + + void CalculateOnIntegrationPoints( + const Variable>& rVariable, + std::vector>& rOutput, + const ProcessInfo& rCurrentProcessInfo) override; + + void CalculateOnIntegrationPoints( + const Variable>& rVariable, + std::vector>& rOutput, + const ProcessInfo& rCurrentProcessInfo) override; + + void CalculateOnIntegrationPoints( + const Variable>& rVariable, + std::vector>& rOutput, + const ProcessInfo& rCurrentProcessInfo) override; + + void CalculateOnIntegrationPoints( + const Variable& rVariable, + std::vector& rOutput, + const ProcessInfo& rCurrentProcessInfo) override; + + void CalculateOnIntegrationPoints( + const Variable& rVariable, + std::vector& rOutput, + const ProcessInfo& rCurrentProcessInfo) override; int Check(const ProcessInfo& rCurrentProcessInfo) const override; diff --git a/applications/SystemIdentificationApplication/custom_sensors/strain_sensor.cpp b/applications/SystemIdentificationApplication/custom_sensors/strain_sensor.cpp index 47161b7066a2..6049236cfe2b 100644 --- a/applications/SystemIdentificationApplication/custom_sensors/strain_sensor.cpp +++ b/applications/SystemIdentificationApplication/custom_sensors/strain_sensor.cpp @@ -41,6 +41,29 @@ StrainSensor::StrainSensor( << "The point " << this->GetLocation() << " is not inside or on the boundary of the geometry of element with id " << mElementId << "."; + // if the element is of type 2d + if (rElement.GetGeometry().WorkingSpaceDimension() == 2) { + switch (mStrainType) { + case StrainType::STRAIN_XX: + mStrainIndex = 0; + break; + case StrainType::STRAIN_YY: + mStrainIndex = 3; + break; + case StrainType::STRAIN_XY: + mStrainIndex = 1; + break; + default: + KRATOS_ERROR << "The element with id = " << rElement.Id() << " is of 2d type, hence only xx, yy, xy strains are allowed."; + } + } else if (rElement.GetGeometry().WorkingSpaceDimension() == 3) { + mStrainIndex = mStrainType; + } else { + KRATOS_ERROR << "Unsupported working space dimension = " + << rElement.GetGeometry().WorkingSpaceDimension() + << " in element with id = " << rElement.Id() << "."; + } + this->SetValue(SENSOR_ELEMENT_ID, static_cast(mElementId)); } @@ -110,7 +133,9 @@ double StrainSensor::CalculateValue(ModelPart& rModelPart) r_element.CalculateOnIntegrationPoints(mrStrainVariable, strains, rModelPart.GetProcessInfo()); for (const auto& strain : strains) { - directional_strain += *(strain.data().begin() + mStrainType); + KRATOS_ERROR_IF(strain.data().size() <= mStrainIndex) + << "The size of the strain " << strain.data().size() << " does not contain the index = " << mStrainIndex << "."; + directional_strain += *(strain.data().begin() + mStrainIndex); } directional_strain /= strains.size(); @@ -138,47 +163,49 @@ void StrainSensor::CalculateGradient( std::vector perturbed_strains, ref_strains; r_element.CalculateOnIntegrationPoints(mrStrainVariable, ref_strains, rProcessInfo); + Element::DofsVectorType elemental_dofs; + r_element.GetDofList(elemental_dofs, rProcessInfo); + const double delta = rProcessInfo[PERTURBATION_SIZE]; auto& r_geometry = r_element.GetGeometry(); + // now only keep the dofs of the first node since we only need the variable + // type to be checked. + const IndexType block_size = elemental_dofs.size() / r_geometry.size(); + elemental_dofs.erase(elemental_dofs.begin() + block_size, elemental_dofs.end()); + IndexType local_index = 0; for (IndexType i_node = 0; i_node < r_geometry.size(); ++i_node) { auto& r_node = r_geometry[i_node]; - if (r_node.HasDofFor(ADJOINT_DISPLACEMENT_X)) { - rResponseGradient[local_index++] = this->CalculateStrainDirectionalSensitivity( - delta, DISPLACEMENT_X, r_node, r_element, - perturbed_strains, ref_strains, rProcessInfo); - } - - if (r_node.HasDofFor(ADJOINT_DISPLACEMENT_Y)) { - rResponseGradient[local_index++] = this->CalculateStrainDirectionalSensitivity( - delta, DISPLACEMENT_Y, r_node, r_element, - perturbed_strains, ref_strains, rProcessInfo); - } - - if (r_node.HasDofFor(ADJOINT_DISPLACEMENT_Z)) { - rResponseGradient[local_index++] = this->CalculateStrainDirectionalSensitivity( - delta, DISPLACEMENT_Z, r_node, r_element, - perturbed_strains, ref_strains, rProcessInfo); - } - - if (r_node.HasDofFor(ADJOINT_ROTATION_X)) { - rResponseGradient[local_index++] = this->CalculateStrainDirectionalSensitivity( - delta, ROTATION_X, r_node, r_element, - perturbed_strains, ref_strains, rProcessInfo); - } - - if (r_node.HasDofFor(ADJOINT_ROTATION_Y)) { - rResponseGradient[local_index++] = this->CalculateStrainDirectionalSensitivity( - delta, ROTATION_Y, r_node, r_element, - perturbed_strains, ref_strains, rProcessInfo); - } - - if (r_node.HasDofFor(ADJOINT_ROTATION_Z)) { - rResponseGradient[local_index++] = this->CalculateStrainDirectionalSensitivity( - delta, ROTATION_Z, r_node, r_element, - perturbed_strains, ref_strains, rProcessInfo); + for (const auto& p_dof : elemental_dofs) { + if (p_dof->GetVariable() == ADJOINT_DISPLACEMENT_X) { + rResponseGradient[local_index++] = this->CalculateStrainDirectionalSensitivity( + delta, DISPLACEMENT_X, r_node, r_element, + perturbed_strains, ref_strains, rProcessInfo); + } else if (p_dof->GetVariable() == ADJOINT_DISPLACEMENT_Y) { + rResponseGradient[local_index++] = this->CalculateStrainDirectionalSensitivity( + delta, DISPLACEMENT_Y, r_node, r_element, + perturbed_strains, ref_strains, rProcessInfo); + } else if (p_dof->GetVariable() == ADJOINT_DISPLACEMENT_Z) { + rResponseGradient[local_index++] = this->CalculateStrainDirectionalSensitivity( + delta, DISPLACEMENT_Z, r_node, r_element, + perturbed_strains, ref_strains, rProcessInfo); + } else if (p_dof->GetVariable() == ADJOINT_ROTATION_X) { + rResponseGradient[local_index++] = this->CalculateStrainDirectionalSensitivity( + delta, ROTATION_X, r_node, r_element, + perturbed_strains, ref_strains, rProcessInfo); + } else if (p_dof->GetVariable() == ADJOINT_ROTATION_Y) { + rResponseGradient[local_index++] = this->CalculateStrainDirectionalSensitivity( + delta, ROTATION_Y, r_node, r_element, + perturbed_strains, ref_strains, rProcessInfo); + } else if (p_dof->GetVariable() == ADJOINT_ROTATION_Z) { + rResponseGradient[local_index++] = this->CalculateStrainDirectionalSensitivity( + delta, ROTATION_Z, r_node, r_element, + perturbed_strains, ref_strains, rProcessInfo); + } else { + KRATOS_ERROR << "Unsupported dof " << p_dof->GetVariable().Name(); + } } } } @@ -285,29 +312,28 @@ void StrainSensor::PrintInfo(std::ostream& rOStream) const void StrainSensor::PrintData(std::ostream& rOStream) const { - PrintInfo(rOStream); rOStream << " Location: " << this->GetLocation() << std::endl; rOStream << " Value: " << this->GetSensorValue() << std::endl; rOStream << " Weight: " << this->GetWeight() << std::endl; rOStream << " Element Id: " << mElementId << std::endl; switch (mStrainType) { case StrainType::STRAIN_XX: - rOStream << " Direction: STRAIN_XX"; + rOStream << " Direction: STRAIN_XX" << std::endl; break; case StrainType::STRAIN_YY: - rOStream << " Direction: STRAIN_YY"; + rOStream << " Direction: STRAIN_YY" << std::endl; break; case StrainType::STRAIN_ZZ: - rOStream << " Direction: STRAIN_ZZ"; + rOStream << " Direction: STRAIN_ZZ" << std::endl; break; case StrainType::STRAIN_XY: - rOStream << " Direction: STRAIN_XY"; + rOStream << " Direction: STRAIN_XY" << std::endl; break; case StrainType::STRAIN_XZ: - rOStream << " Direction: STRAIN_XZ"; + rOStream << " Direction: STRAIN_XZ" << std::endl; break; case StrainType::STRAIN_YZ: - rOStream << " Direction: STRAIN_YZ"; + rOStream << " Direction: STRAIN_YZ" << std::endl; break; } DataValueContainer::PrintData(rOStream); @@ -341,8 +367,8 @@ double StrainSensor::CalculateStrainDirectionalSensitivity( double strain_sensitivity = 0.0; for (IndexType j = 0; j < rPerturbedStrains.size(); ++j) { - strain_sensitivity += (*(rPerturbedStrains[j].data().begin() + mStrainType) - - *(rRefStrains[j].data().begin() + mStrainType)) / + strain_sensitivity += (*(rPerturbedStrains[j].data().begin() + mStrainIndex) - + *(rRefStrains[j].data().begin() + mStrainIndex)) / Perturbation; } diff --git a/applications/SystemIdentificationApplication/custom_sensors/strain_sensor.h b/applications/SystemIdentificationApplication/custom_sensors/strain_sensor.h index b49e8654dacc..68c9536bd3ef 100644 --- a/applications/SystemIdentificationApplication/custom_sensors/strain_sensor.h +++ b/applications/SystemIdentificationApplication/custom_sensors/strain_sensor.h @@ -167,6 +167,8 @@ class KRATOS_API(SYSTEM_IDENTIFICATION_APPLICATION) StrainSensor : public Sensor StrainType mStrainType; + IndexType mStrainIndex; + Point mLocalPoint; const Variable& mrStrainVariable; diff --git a/applications/SystemIdentificationApplication/tests/auxiliary_files/solid_material_properties.json b/applications/SystemIdentificationApplication/tests/auxiliary_files/solid_material_properties.json new file mode 100644 index 000000000000..a932c3d35b42 --- /dev/null +++ b/applications/SystemIdentificationApplication/tests/auxiliary_files/solid_material_properties.json @@ -0,0 +1,19 @@ +{ + "properties": [ + { + "model_part_name": "Test", + "properties_id": 2, + "Material": { + "constitutive_law": { + "name": "LinearElastic3DLaw" + }, + "Variables": { + "DENSITY": 1.0, + "YOUNG_MODULUS": 30000000000.0, + "POISSON_RATIO": 0.29 + }, + "Tables": {} + } + } + ] +} \ No newline at end of file diff --git a/applications/SystemIdentificationApplication/tests/responses/auxiliary_files_2/StructuralMaterials.json b/applications/SystemIdentificationApplication/tests/responses/auxiliary_files_2/StructuralMaterials.json new file mode 100644 index 000000000000..edca5cfa2620 --- /dev/null +++ b/applications/SystemIdentificationApplication/tests/responses/auxiliary_files_2/StructuralMaterials.json @@ -0,0 +1,18 @@ +{ + "properties" : [{ + "model_part_name" : "Structure.Parts_Solid_full", + "properties_id" : 1, + "Material" : { + "constitutive_law" : { + "name" : "LinearElastic3DLaw" + }, + "Variables" : { + "THICKNESS" : 0.02, + "DENSITY" : 7850.0, + "YOUNG_MODULUS" : 20.0, + "POISSON_RATIO" : 0.29 + }, + "Tables" : null + } + }] +} diff --git a/applications/SystemIdentificationApplication/tests/responses/auxiliary_files_2/Structure.mdpa b/applications/SystemIdentificationApplication/tests/responses/auxiliary_files_2/Structure.mdpa new file mode 100644 index 000000000000..53c1a40e0abe --- /dev/null +++ b/applications/SystemIdentificationApplication/tests/responses/auxiliary_files_2/Structure.mdpa @@ -0,0 +1,79 @@ +Begin ModelPartData +// VARIABLE_NAME value +End ModelPartData + +Begin Properties 0 +End Properties +Begin Nodes + 1 0.0000000000 0.0000000000 0.0000000000 + 2 1.0000000000 0.0000000000 0.0000000000 + 3 1.0000000000 1.0000000000 0.0000000000 + 4 0.0000000000 1.0000000000 0.0000000000 + 5 0.5000000000 0.5000000000 1.0000000000 +End Nodes + + +Begin Elements SmallDisplacementElement3D4N// GUI group identifier: volume + 1 0 1 2 3 5 + 2 0 1 3 4 5 +End Elements + +Begin Conditions SurfaceLoadCondition3D3N// GUI group identifier: right + 2 0 2 3 5 +End Conditions + +Begin SubModelPart Parts_Mat_1 + Begin SubModelPartNodes + 1 + 2 + 3 + 5 + End SubModelPartNodes + Begin SubModelPartElements + 2 + End SubModelPartElements +End SubModelPart + +Begin SubModelPart Parts_Solid_full + Begin SubModelPartNodes + 1 + 2 + 3 + 4 + 5 + End SubModelPartNodes + Begin SubModelPartElements + 1 + 2 + End SubModelPartElements +End SubModelPart + +Begin SubModelPart Parts_Solid_volume // Group volume // Subtree Parts_Solid + Begin SubModelPartNodes + 1 + 2 + 3 + 5 + End SubModelPartNodes + Begin SubModelPartElements + 1 + End SubModelPartElements +End SubModelPart +Begin SubModelPart DISPLACEMENT_left // Group left // Subtree DISPLACEMENT + Begin SubModelPartNodes + 1 + 2 + 3 + 4 + End SubModelPartNodes +End SubModelPart +Begin SubModelPart LineLoad2D_right // Group right // Subtree LineLoad2D + Begin SubModelPartNodes + 2 + 3 + 5 + End SubModelPartNodes + Begin SubModelPartConditions + 2 + End SubModelPartConditions +End SubModelPart diff --git a/applications/SystemIdentificationApplication/tests/responses/auxiliary_files_2/adjoint_material_parameters.json b/applications/SystemIdentificationApplication/tests/responses/auxiliary_files_2/adjoint_material_parameters.json new file mode 100644 index 000000000000..85b9c6d34737 --- /dev/null +++ b/applications/SystemIdentificationApplication/tests/responses/auxiliary_files_2/adjoint_material_parameters.json @@ -0,0 +1,3 @@ +{ + "properties": [] +} \ No newline at end of file diff --git a/applications/SystemIdentificationApplication/tests/responses/auxiliary_files_2/adjoint_project_parameters_p_norm.json b/applications/SystemIdentificationApplication/tests/responses/auxiliary_files_2/adjoint_project_parameters_p_norm.json new file mode 100644 index 000000000000..17de29a526c1 --- /dev/null +++ b/applications/SystemIdentificationApplication/tests/responses/auxiliary_files_2/adjoint_project_parameters_p_norm.json @@ -0,0 +1,77 @@ +{ + "problem_data": { + "problem_name": "benchmark_2", + "parallel_type": "OpenMP", + "start_time": 0, + "end_time": 1, + "echo_level": 0 + }, + "sensor_settings": { + "perturbation_size": 1e-8, + "adapt_perturbation_size": true, + "p_coefficient": 4, + "@include_json": "auxiliary_files_2/sensor_data.json" + }, + "solver_settings": { + "solver_type": "adjoint_static", + "analysis_type": "linear", + "model_part_name": "AdjointStructure", + "domain_size": 3, + "time_stepping": { + "time_step": 1.0 + }, + "compute_reactions": false, + "move_mesh_flag": false, + "sensitivity_settings": { + "sensitivity_model_part_name": "Parts_Solid_full", + "element_data_value_sensitivity_variables": [ + "YOUNG_MODULUS" + ], + "build_mode": "static" + }, + "echo_level": 0, + "rotation_dofs": true, + "model_import_settings": { + "input_type": "use_input_model_part" + }, + "material_import_settings": { + "materials_filename": "auxiliary_files_2/adjoint_material_parameters.json" + }, + "linear_solver_settings": { + "solver_type": "amgcl" + } + }, + "processes": { + "constraints_process_list": [ + { + "python_module": "assign_vector_variable_process", + "kratos_module": "KratosMultiphysics", + "help": "This process fixes the selected components of a given vector variable", + "process_name": "AssignVectorVariableProcess", + "Parameters": { + "mesh_id": 0, + "model_part_name": "AdjointStructure.DISPLACEMENT_left", + "variable_name": "ADJOINT_DISPLACEMENT", + "value": [ + 0.0, + 0.0, + 0.0 + ], + "constrained": [ + true, + true, + true + ], + "interval": [ + 0.0, + "End" + ] + } + } + ], + "loads_process_list": [], + "list_other_processes": [] + }, + "output_processes": {}, + "analysis_stage": "KratosMultiphysics.SystemIdentificationApplication.sensor_sensitivity_solvers.sensor_sensitivity_static_analysis" +} \ No newline at end of file diff --git a/applications/SystemIdentificationApplication/tests/responses/auxiliary_files_2/damaged_system/MainKratos.py b/applications/SystemIdentificationApplication/tests/responses/auxiliary_files_2/damaged_system/MainKratos.py new file mode 100644 index 000000000000..bb644f15e8a0 --- /dev/null +++ b/applications/SystemIdentificationApplication/tests/responses/auxiliary_files_2/damaged_system/MainKratos.py @@ -0,0 +1,16 @@ +import KratosMultiphysics as Kratos +from KratosMultiphysics.StructuralMechanicsApplication.structural_mechanics_analysis import StructuralMechanicsAnalysis + +class CustomStructuralMechanicsAnalysis(StructuralMechanicsAnalysis): + def FinalizeSolutionStep(self): + for element in self._GetSolver().GetComputingModelPart().Elements: + element.SetValue(Kratos.YOUNG_MODULUS, element.Properties[Kratos.YOUNG_MODULUS]) + +if __name__ == "__main__": + model = Kratos.Model() + + with open("primal_project_parameters.json", "r") as file_input: + parameters = Kratos.Parameters(file_input.read()) + + analysis = CustomStructuralMechanicsAnalysis(model, parameters) + analysis.Run() \ No newline at end of file diff --git a/applications/SystemIdentificationApplication/tests/responses/auxiliary_files_2/damaged_system/StructuralMaterials.json b/applications/SystemIdentificationApplication/tests/responses/auxiliary_files_2/damaged_system/StructuralMaterials.json new file mode 100644 index 000000000000..0ca7036968b7 --- /dev/null +++ b/applications/SystemIdentificationApplication/tests/responses/auxiliary_files_2/damaged_system/StructuralMaterials.json @@ -0,0 +1,34 @@ +{ + "properties" : [ + { + "model_part_name" : "Structure.Parts_Solid_full", + "properties_id" : 1, + "Material" : { + "constitutive_law" : { + "name" : "LinearElastic3DLaw" + }, + "Variables" : { + "DENSITY" : 7850.0, + "YOUNG_MODULUS" : 20.0, + "POISSON_RATIO" : 0.29 + }, + "Tables" : null + } + }, + { + "model_part_name" : "Structure.Parts_Mat_1", + "properties_id" : 2, + "Material" : { + "constitutive_law" : { + "name" : "LinearElastic3DLaw" + }, + "Variables" : { + "DENSITY" : 7850.0, + "YOUNG_MODULUS" : 5.0, + "POISSON_RATIO" : 0.29 + }, + "Tables" : null + } + } +] +} diff --git a/applications/SystemIdentificationApplication/tests/responses/auxiliary_files_2/damaged_system/primal_project_parameters.json b/applications/SystemIdentificationApplication/tests/responses/auxiliary_files_2/damaged_system/primal_project_parameters.json new file mode 100644 index 000000000000..23ac76c4a98f --- /dev/null +++ b/applications/SystemIdentificationApplication/tests/responses/auxiliary_files_2/damaged_system/primal_project_parameters.json @@ -0,0 +1,137 @@ +{ + "problem_data": { + "problem_name": "gg", + "parallel_type": "OpenMP", + "echo_level": 1, + "start_time": 0.0, + "end_time": 1.0 + }, + "solver_settings": { + "time_stepping": { + "time_step": 1.1 + }, + "solver_type": "Static", + "model_part_name": "Structure", + "domain_size": 2, + "echo_level": 0, + "analysis_type": "non_linear", + "model_import_settings": { + "input_type": "mdpa", + "input_filename": "../Structure" + }, + "material_import_settings": { + "materials_filename": "StructuralMaterials.json" + }, + "line_search": false, + "convergence_criterion": "residual_criterion", + "displacement_relative_tolerance": 0.0001, + "displacement_absolute_tolerance": 1e-9, + "residual_relative_tolerance": 0.0001, + "residual_absolute_tolerance": 1e-9, + "max_iteration": 10, + "use_old_stiffness_in_first_iteration": false, + "rotation_dofs": false, + "volumetric_strain_dofs": false, + "move_mesh_flag": false + }, + "processes": { + "constraints_process_list": [ + { + "python_module": "assign_vector_variable_process", + "kratos_module": "KratosMultiphysics", + "help": "This process fixes the selected components of a given vector variable", + "process_name": "AssignVectorVariableProcess", + "Parameters": { + "mesh_id": 0, + "model_part_name": "Structure.DISPLACEMENT_left", + "variable_name": "DISPLACEMENT", + "value": [ + 0.0, + 0.0, + 0.0 + ], + "constrained": [ + true, + true, + true + ], + "interval": [ + 0.0, + "End" + ] + } + } + ], + "loads_process_list": [ + { + "python_module": "assign_vector_by_direction_to_condition_process", + "kratos_module": "KratosMultiphysics", + "help": "This process sets a vector variable value over a condition", + "check": "DirectorVectorNonZero direction", + "process_name": "AssignModulusAndDirectionToConditionsProcess", + "Parameters": { + "model_part_name": "Structure.LineLoad2D_right", + "variable_name": "SURFACE_LOAD", + "modulus": 0.01, + "direction": [ + 1.0, + 0.0, + 0.0 + ], + "interval": [ + 0.0, + "End" + ] + } + } + ], + "list_other_processes": [] + }, + "output_processes": { + "vtk_output": [ + { + "Parameters": { + "condition_data_value_variables": [], + "element_data_value_variables": [ + "YOUNG_MODULUS" + ], + "file_format": "binary", + "model_part_name": "Structure", + "nodal_data_value_variables": [], + "nodal_solution_step_data_variables": [ + "DISPLACEMENT", + "REACTION", + "ROTATION" + ], + "gauss_point_variables_extrapolated_to_nodes": [ + "VON_MISES_STRESS", + "GREEN_LAGRANGE_STRAIN_VECTOR" + ], + "output_control_type": "step", + "output_interval": 1, + "output_path": "Output", + "output_precision": 7, + "output_sub_model_parts": true, + "save_output_files_in_folder": true + }, + "help": "This process writes postprocessing files for Paraview", + "kratos_module": "KratosMultiphysics", + "process_name": "VtkOutputProcess", + "python_module": "vtk_output_process" + } + ], + "sensor_output": [ + { + "Parameters": { + "model_part_name": "Structure", + "output_file_name": "../measured_data.csv", + "@include_json": "../sensor_data.json" + }, + "kratos_module": "KratosMultiphysics.SystemIdentificationApplication.processes", + "process_name": "SensorOutputProcess", + "python_module": "sensor_output_process" + } + ] + }, + "analysis_stage": "KratosMultiphysics.StructuralMechanicsApplication.structural_mechanics_analysis" +} diff --git a/applications/SystemIdentificationApplication/tests/responses/auxiliary_files_2/measured_data.csv b/applications/SystemIdentificationApplication/tests/responses/auxiliary_files_2/measured_data.csv new file mode 100644 index 000000000000..0c9b1120cc0f --- /dev/null +++ b/applications/SystemIdentificationApplication/tests/responses/auxiliary_files_2/measured_data.csv @@ -0,0 +1,3 @@ +#,type,name,location_0,location_1,location_2,value +1,strain_sensor,strain_xx_1,0.2,0.2,0.2,0.0005769055381949457 +2,strain_sensor,strain_xx_2,0.8,0.8,0.2,0.0005769055381949457 diff --git a/applications/SystemIdentificationApplication/tests/responses/auxiliary_files_2/optimization_parameters_p_norm.json b/applications/SystemIdentificationApplication/tests/responses/auxiliary_files_2/optimization_parameters_p_norm.json new file mode 100644 index 000000000000..a82b33c4745f --- /dev/null +++ b/applications/SystemIdentificationApplication/tests/responses/auxiliary_files_2/optimization_parameters_p_norm.json @@ -0,0 +1,133 @@ +{ + "problem_data": { + "parallel_type": "OpenMP", + "echo_level": 0 + }, + "model_parts": [ + { + "type": "mdpa_model_part_controller", + "settings": { + "model_part_name": "AdjointStructure", + "domain_size": 3, + "input_filename": "auxiliary_files_2/Structure" + } + }, + { + "type": "connectivity_preserving_model_part_controller", + "settings": { + "transformation_settings": [ + { + "source_model_part_name": "AdjointStructure", + "destination_model_part_name": "Structure", + "destination_element_name": "SmallDisplacementElement3D4N", + "destination_condition_name": "SurfaceLoadCondition3D3N" + } + ] + } + } + ], + "analyses": [ + { + "name": "Structure_static", + "type": "kratos_analysis_execution_policy", + "settings": { + "model_part_names": [ + "Structure" + ], + "analysis_module": "KratosMultiphysics.StructuralMechanicsApplication", + "analysis_type": "StructuralMechanicsAnalysis", + "analysis_settings": { + "@include_json": "auxiliary_files_2/primal_project_parameters.json" + } + } + } + ], + "responses": [ + { + "name": "damage_response", + "type": "damage_detection_response", + "module": "KratosMultiphysics.SystemIdentificationApplication.responses", + "settings": { + "evaluated_model_part_names": [ + "AdjointStructure" + ], + "adjoint_parameters": { + "@include_json": "auxiliary_files_2/adjoint_project_parameters_p_norm.json" + }, + "test_analysis_list": [ + { + "primal_analysis_name": "Structure_static", + "sensor_measurement_csv_file": "auxiliary_files_2/measured_data.csv", + "weight": 1.0 + } + ] + } + } + ], + "controls": [ + { + "name": "material_control", + "type": "material_properties_control", + "module": "KratosMultiphysics.SystemIdentificationApplication.controls", + "settings": { + "model_part_names": [ + { + "primal_model_part_name": "Structure", + "adjoint_model_part_name": "AdjointStructure" + } + ], + "control_variable_name": "YOUNG_MODULUS", + "control_variable_bounds": [ + 0.0, + 30000000000.0 + ], + "filter_settings": { + "filter_type": "explicit_filter", + "filter_radius_settings": { + "filter_radius": 5.0, + "filter_radius_type": "constant" + } + } + } + } + ], + "algorithm_settings": { + "type": "algorithm_steepest_descent", + "settings": { + "echo_level": 0, + "line_search": { + "type": "const_step", + "init_step": 0.01, + "gradient_scaling": "inf_norm" + }, + "conv_settings": { + "type": "max_iter", + "max_iter": 5 + } + }, + "controls": [ + "material_control" + ], + "objective": { + "response_name": "damage_response", + "type": "minimization", + "scaling": 1.0 + } + }, + "processes": { + "kratos_processes": {}, + "optimization_data_processes": { + "output_processes": [ + { + "type": "optimization_problem_ascii_output_process", + "module": "KratosMultiphysics.OptimizationApplication.processes", + "settings": { + "output_file_name": "auxiliary_files_2/summary.csv", + "write_kratos_version": false, + "write_time_stamp": false + } + } + ] + } + } +} \ No newline at end of file diff --git a/applications/SystemIdentificationApplication/tests/responses/auxiliary_files_2/primal_project_parameters.json b/applications/SystemIdentificationApplication/tests/responses/auxiliary_files_2/primal_project_parameters.json new file mode 100644 index 000000000000..0ac57c2c286e --- /dev/null +++ b/applications/SystemIdentificationApplication/tests/responses/auxiliary_files_2/primal_project_parameters.json @@ -0,0 +1,92 @@ +{ + "problem_data": { + "problem_name": "gg", + "parallel_type": "OpenMP", + "echo_level": 1, + "start_time": 0.0, + "end_time": 1.0 + }, + "solver_settings": { + "time_stepping": { + "time_step": 1.1 + }, + "solver_type": "Static", + "model_part_name": "Structure", + "domain_size": 2, + "echo_level": 0, + "analysis_type": "non_linear", + "model_import_settings": { + "input_type": "use_input_model_part" + }, + "material_import_settings": { + "materials_filename": "auxiliary_files_2/StructuralMaterials.json" + }, + "line_search": false, + "convergence_criterion": "residual_criterion", + "displacement_relative_tolerance": 0.0001, + "displacement_absolute_tolerance": 1e-9, + "residual_relative_tolerance": 0.0001, + "residual_absolute_tolerance": 1e-9, + "max_iteration": 10, + "use_old_stiffness_in_first_iteration": false, + "rotation_dofs": true, + "volumetric_strain_dofs": false, + "move_mesh_flag": false + }, + "processes": { + "constraints_process_list": [ + { + "python_module": "assign_vector_variable_process", + "kratos_module": "KratosMultiphysics", + "help": "This process fixes the selected components of a given vector variable", + "process_name": "AssignVectorVariableProcess", + "Parameters": { + "mesh_id": 0, + "model_part_name": "Structure.DISPLACEMENT_left", + "variable_name": "DISPLACEMENT", + "value": [ + 0.0, + 0.0, + 0.0 + ], + "constrained": [ + true, + true, + true + ], + "interval": [ + 0.0, + "End" + ] + } + } + ], + "loads_process_list": [ + { + "python_module": "assign_vector_by_direction_to_condition_process", + "kratos_module": "KratosMultiphysics", + "help": "This process sets a vector variable value over a condition", + "check": "DirectorVectorNonZero direction", + "process_name": "AssignModulusAndDirectionToConditionsProcess", + "Parameters": { + "model_part_name": "Structure.LineLoad2D_right", + "variable_name": "LINE_LOAD", + "modulus": 0.01, + "direction": [ + 1.0, + 0.0, + 0.0 + ], + "interval": [ + 0.0, + "End" + ] + } + } + ], + "list_other_processes": [] + }, + "output_processes": { + }, + "analysis_stage": "KratosMultiphysics.StructuralMechanicsApplication.structural_mechanics_analysis" +} diff --git a/applications/SystemIdentificationApplication/tests/responses/auxiliary_files_2/sensor_data.json b/applications/SystemIdentificationApplication/tests/responses/auxiliary_files_2/sensor_data.json new file mode 100644 index 000000000000..7812443114d8 --- /dev/null +++ b/applications/SystemIdentificationApplication/tests/responses/auxiliary_files_2/sensor_data.json @@ -0,0 +1,30 @@ +{ + "list_of_sensors": [ + { + "location": [ + 0.2, + 0.2, + 0.2 + ], + "name": "strain_xx_1", + "strain_type": "strain_xz", + "strain_variable": "GREEN_LAGRANGE_STRAIN_TENSOR", + "type": "strain_sensor", + "variable_data": {}, + "weight": 1.0 + }, + { + "location": [ + 0.8, + 0.8, + 0.2 + ], + "name": "strain_xx_2", + "strain_type": "strain_xz", + "strain_variable": "GREEN_LAGRANGE_STRAIN_TENSOR", + "type": "strain_sensor", + "variable_data": {}, + "weight": 1.0 + } + ] +} \ No newline at end of file diff --git a/applications/SystemIdentificationApplication/tests/responses/test_damage_response.py b/applications/SystemIdentificationApplication/tests/responses/test_damage_response.py index 020ce757fe68..aecdfe769fd5 100644 --- a/applications/SystemIdentificationApplication/tests/responses/test_damage_response.py +++ b/applications/SystemIdentificationApplication/tests/responses/test_damage_response.py @@ -152,5 +152,36 @@ def test_DamageResponse(self): self.assertAlmostEqual(gradients[index], sensitivity, 6) element.Properties[Kratos.YOUNG_MODULUS] -= delta +class TestDamageDetectionResponseStrainSensor(kratos_unittest.TestCase): + def test_DamageResponse(self): + with kratos_unittest.WorkFolderScope(".", __file__): + with open("auxiliary_files_2/optimization_parameters_p_norm.json", "r") as file_input: + parameters = Kratos.Parameters(file_input.read()) + + model = Kratos.Model() + analysis = OptimizationAnalysis(model, parameters) + + analysis.Initialize() + analysis.Check() + objective: ResponseRoutine = analysis.optimization_problem.GetComponent("damage_response", ResponseRoutine) + + var = objective.GetRequiredPhysicalGradients() + response = objective.GetReponse() + model_part = response.GetInfluencingModelPart() + + ref_value = response.CalculateValue() + self.assertAlmostEqual(ref_value, 1.9789595600760277e-07, 6) + + response.CalculateGradient(var) + + gradients = var[Kratos.YOUNG_MODULUS].Evaluate() + + delta = 1e-8 + for index, element in enumerate(model_part.Elements): + element.Properties[Kratos.YOUNG_MODULUS] += delta + sensitivity = ((response.CalculateValue() - ref_value) / delta) + self.assertAlmostEqual(gradients[index], sensitivity, 6) + element.Properties[Kratos.YOUNG_MODULUS] -= delta + if __name__ == "__main__": kratos_unittest.main() \ No newline at end of file diff --git a/applications/SystemIdentificationApplication/tests/test_SystemIdentificationApplication.py b/applications/SystemIdentificationApplication/tests/test_SystemIdentificationApplication.py index ad5099b1d248..fea9d0ba84af 100644 --- a/applications/SystemIdentificationApplication/tests/test_SystemIdentificationApplication.py +++ b/applications/SystemIdentificationApplication/tests/test_SystemIdentificationApplication.py @@ -14,7 +14,8 @@ def AssembleTestSuites(): smallSuite = suites['small'] smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_adjoint_sensors.TestDisplacementSensor])) - smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_adjoint_sensors.TestStrainSensor])) + smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_adjoint_sensors.TestStrainSensorShell])) + smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_adjoint_sensors.TestStrainSensorSolids])) smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_sensor_output_process.TestSensorOutputProcess])) smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_smooth_clamper.TestSmoothClamper])) smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_mask_utils.TestMaskUtils])) diff --git a/applications/SystemIdentificationApplication/tests/test_adjoint_sensors.py b/applications/SystemIdentificationApplication/tests/test_adjoint_sensors.py index 1ad9b753c305..bc397113656f 100644 --- a/applications/SystemIdentificationApplication/tests/test_adjoint_sensors.py +++ b/applications/SystemIdentificationApplication/tests/test_adjoint_sensors.py @@ -1,5 +1,6 @@ from math import sqrt import KratosMultiphysics as Kratos +import KratosMultiphysics.OptimizationApplication as KratosOA import KratosMultiphysics.SystemIdentificationApplication as KratosSI import KratosMultiphysics.StructuralMechanicsApplication as KratosStruct import KratosMultiphysics.KratosUnittest as UnitTest @@ -191,31 +192,39 @@ def test_CalculateGradient(self): element: Kratos.Element = self.model_part.GetElement(sensor.GetValue(KratosSI.SENSOR_ELEMENT_ID)) sensor.CalculateGradient(element, residual_matrix, response_sensitivities, self.model_part.ProcessInfo) - for i, node in enumerate(element.GetGeometry()): + for j, node in enumerate(element.GetGeometry()): node.SetSolutionStepValue(Kratos.DISPLACEMENT_X, node.GetSolutionStepValue(Kratos.DISPLACEMENT_X) + delta) perturbed_value = sensor.CalculateValue(self.model_part) sensitivity = (perturbed_value - ref_value) / delta - self.assertAlmostEqual(sensitivity, response_sensitivities[i * 6]) + self.assertAlmostEqual(sensitivity, response_sensitivities[j * 6]) node.SetSolutionStepValue(Kratos.DISPLACEMENT_X, node.GetSolutionStepValue(Kratos.DISPLACEMENT_X) - delta) node.SetSolutionStepValue(Kratos.DISPLACEMENT_Y, node.GetSolutionStepValue(Kratos.DISPLACEMENT_Y) + delta) perturbed_value = sensor.CalculateValue(self.model_part) sensitivity = (perturbed_value - ref_value) / delta - self.assertAlmostEqual(sensitivity, response_sensitivities[i * 6 + 1]) + self.assertAlmostEqual(sensitivity, response_sensitivities[j * 6 + 1]) node.SetSolutionStepValue(Kratos.DISPLACEMENT_Y, node.GetSolutionStepValue(Kratos.DISPLACEMENT_Y) - delta) -class TestStrainSensor(UnitTest.TestCase): +class TestStrainSensorShell(UnitTest.TestCase): @classmethod def setUpClass(cls) -> None: cls.model = Kratos.Model() cls.model_part = cls.model.CreateModelPart("Test") + cls.adjoint_model_part = cls.model.CreateModelPart("TestAdjoint") + cls.model_part.ProcessInfo[KratosStruct.PERTURBATION_SIZE] = 1e-10 + cls.model_part.AddNodalSolutionStepVariable(Kratos.DISPLACEMENT) - cls.model_part.AddNodalSolutionStepVariable(KratosStruct.ADJOINT_DISPLACEMENT) cls.model_part.AddNodalSolutionStepVariable(Kratos.ROTATION) + cls.model_part.AddNodalSolutionStepVariable(KratosStruct.ADJOINT_DISPLACEMENT) cls.model_part.AddNodalSolutionStepVariable(KratosStruct.ADJOINT_ROTATION) + cls.adjoint_model_part.AddNodalSolutionStepVariable(Kratos.DISPLACEMENT) + cls.adjoint_model_part.AddNodalSolutionStepVariable(Kratos.ROTATION) + cls.adjoint_model_part.AddNodalSolutionStepVariable(KratosStruct.ADJOINT_DISPLACEMENT) + cls.adjoint_model_part.AddNodalSolutionStepVariable(KratosStruct.ADJOINT_ROTATION) + cls.model_part.CreateNewNode(1, 0.0, 0.0, 0.0) cls.model_part.CreateNewNode(2, 1.0, 0.0, 0.0) cls.model_part.CreateNewNode(3, 1.0, 1.0, 0.0) @@ -230,9 +239,16 @@ def setUpClass(cls) -> None: material_settings = Kratos.Parameters("""{"Parameters": {"materials_filename": "auxiliary_files/shell_material_properties.json"}}""") Kratos.ReadMaterialsUtility(material_settings, cls.model) + KratosOA.OptimizationUtils.CreateEntitySpecificPropertiesForContainer(cls.model_part, cls.model_part.Elements, False) + + Kratos.ConnectivityPreserveModeler().GenerateModelPart(cls.model_part, cls.adjoint_model_part, "AdjointFiniteDifferencingShellThinElement3D3N") + for element in cls.model_part.Elements: element.Initialize(cls.model_part.ProcessInfo) + for element in cls.adjoint_model_part.Elements: + element.Initialize(cls.adjoint_model_part.ProcessInfo) + Kratos.VariableUtils().AddDof(Kratos.DISPLACEMENT_X, cls.model_part) Kratos.VariableUtils().AddDof(Kratos.DISPLACEMENT_Y, cls.model_part) Kratos.VariableUtils().AddDof(Kratos.DISPLACEMENT_Z, cls.model_part) @@ -310,9 +326,9 @@ def test_CalculateGradient(self): ref_value = self.ref_values[i] delta = 1e-5 - element: Kratos.Element = self.model_part.GetElement(sensor.GetValue(KratosSI.SENSOR_ELEMENT_ID)) - sensor.CalculateGradient(element, residual_matrix, response_sensitivities, self.model_part.ProcessInfo) - for i, node in enumerate(element.GetGeometry()): + adjoint_element: Kratos.Element = self.adjoint_model_part.GetElement(sensor.GetValue(KratosSI.SENSOR_ELEMENT_ID)) + sensor.CalculateGradient(adjoint_element, residual_matrix, response_sensitivities, self.model_part.ProcessInfo) + for i, node in enumerate(adjoint_element.GetGeometry()): node.SetSolutionStepValue(Kratos.DISPLACEMENT_X, node.GetSolutionStepValue(Kratos.DISPLACEMENT_X) + delta) perturbed_value = sensor.CalculateValue(self.model_part) sensitivity = (perturbed_value - ref_value) / delta @@ -326,5 +342,225 @@ def test_CalculateGradient(self): node.SetSolutionStepValue(Kratos.DISPLACEMENT_Y, node.GetSolutionStepValue(Kratos.DISPLACEMENT_Y) - delta) +class TestStrainSensorSolids(UnitTest.TestCase): + @classmethod + def setUpClass(cls) -> None: + cls.model = Kratos.Model() + cls.model_part = cls.model.CreateModelPart("Test") + cls.adjoint_model_part = cls.model.CreateModelPart("TestAdjoint") + + cls.model_part.ProcessInfo[Kratos.DOMAIN_SIZE] = 3 + cls.model_part.ProcessInfo[KratosStruct.PERTURBATION_SIZE] = 1e-10 + cls.model_part.AddNodalSolutionStepVariable(Kratos.DISPLACEMENT) + cls.model_part.AddNodalSolutionStepVariable(KratosStruct.ADJOINT_DISPLACEMENT) + + cls.adjoint_model_part.AddNodalSolutionStepVariable(Kratos.DISPLACEMENT) + cls.adjoint_model_part.AddNodalSolutionStepVariable(KratosStruct.ADJOINT_DISPLACEMENT) + + cls.model_part.CreateNewNode( 1, 0.0, 0.0, 0.0) + cls.model_part.CreateNewNode( 2, 1.0, 0.0, 0.0) + cls.model_part.CreateNewNode( 3, 1.0, 1.0, 0.0) + cls.model_part.CreateNewNode( 4, 0.0, 1.0, 0.0) + cls.model_part.CreateNewNode( 5, 0.0, 0.0, 1.0) + cls.model_part.CreateNewNode( 6, 1.0, 0.0, 1.0) + cls.model_part.CreateNewNode( 7, 1.0, 1.0, 1.0) + cls.model_part.CreateNewNode( 8, 0.0, 1.0, 1.0) + cls.model_part.CreateNewNode( 9, 2.0, 0.0, 0.0) + cls.model_part.CreateNewNode(10, 2.0, 1.0, 0.0) + cls.model_part.CreateNewNode(11, 2.0, 0.0, 1.0) + cls.model_part.CreateNewNode(12, 2.0, 1.0, 1.0) + + prop = cls.model_part.CreateNewProperties(1) + + cls.model_part.CreateNewElement("SmallDisplacementElement3D8N", 1, [1, 2, 3, 4, 5, 6, 7, 8], prop) + cls.model_part.CreateNewElement("SmallDisplacementElement3D8N", 2, [2, 9, 10, 3, 6, 11, 12, 7], prop) + + # set the constitutive laws + material_settings = Kratos.Parameters("""{"Parameters": {"materials_filename": "auxiliary_files/solid_material_properties.json"}}""") + Kratos.ReadMaterialsUtility(material_settings, cls.model) + + KratosOA.OptimizationUtils.CreateEntitySpecificPropertiesForContainer(cls.model_part, cls.model_part.Elements, False) + + Kratos.ConnectivityPreserveModeler().GenerateModelPart(cls.model_part, cls.adjoint_model_part, "AdjointFiniteDifferencingSmallDisplacementElement3D8N") + + for element in cls.model_part.Elements: + element.Initialize(cls.model_part.ProcessInfo) + + for element in cls.adjoint_model_part.Elements: + element.Initialize(cls.adjoint_model_part.ProcessInfo) + + Kratos.VariableUtils().AddDof(Kratos.DISPLACEMENT_X, cls.model_part) + Kratos.VariableUtils().AddDof(Kratos.DISPLACEMENT_Y, cls.model_part) + Kratos.VariableUtils().AddDof(Kratos.DISPLACEMENT_Z, cls.model_part) + Kratos.VariableUtils().AddDof(KratosStruct.ADJOINT_DISPLACEMENT_X, cls.model_part) + Kratos.VariableUtils().AddDof(KratosStruct.ADJOINT_DISPLACEMENT_Y, cls.model_part) + Kratos.VariableUtils().AddDof(KratosStruct.ADJOINT_DISPLACEMENT_Z, cls.model_part) + + for node in cls.model_part.Nodes: + node.SetSolutionStepValue(Kratos.DISPLACEMENT, [node.Id, node.Id + 1, node.Id + 2]) + node.SetValue(Kratos.DISPLACEMENT, node.GetSolutionStepValue(Kratos.DISPLACEMENT)) + + parameters = [ + Kratos.Parameters("""{ + "type" : "strain_sensor", + "name" : "strain_x_1", + "value" : 0, + "location" : [0.5, 0.5, 0.5], + "strain_type" : "strain_xx", + "strain_variable": "GREEN_LAGRANGE_STRAIN_TENSOR", + "weight" : 1.0, + "variable_data" : {} + }"""), + Kratos.Parameters("""{ + "type" : "strain_sensor", + "name" : "strain_x_2", + "value" : 0, + "location" : [1.5, 0.5, 0.5], + "strain_type" : "strain_xx", + "strain_variable": "GREEN_LAGRANGE_STRAIN_TENSOR", + "weight" : 1.0, + "variable_data" : {} + }"""), + Kratos.Parameters("""{ + "type" : "strain_sensor", + "name" : "strain_y_1", + "value" : 0, + "location" : [0.5, 0.5, 0.5], + "strain_type" : "strain_yy", + "strain_variable": "GREEN_LAGRANGE_STRAIN_TENSOR", + "weight" : 1.0, + "variable_data" : {} + }"""), + Kratos.Parameters("""{ + "type" : "strain_sensor", + "name" : "strain_y_2", + "value" : 0, + "location" : [1.5, 0.5, 0.5], + "strain_type" : "strain_yy", + "strain_variable": "GREEN_LAGRANGE_STRAIN_TENSOR", + "weight" : 1.0, + "variable_data" : {} + }"""), + Kratos.Parameters("""{ + "type" : "strain_sensor", + "name" : "strain_xy_1", + "value" : 0, + "location" : [0.5, 0.5, 0.5], + "strain_type" : "strain_xy", + "strain_variable": "GREEN_LAGRANGE_STRAIN_TENSOR", + "weight" : 1.0, + "variable_data" : {} + }"""), + Kratos.Parameters("""{ + "type" : "strain_sensor", + "name" : "strain_xy_2", + "value" : 0, + "location" : [1.5, 0.5, 0.5], + "strain_type" : "strain_xy", + "strain_variable": "GREEN_LAGRANGE_STRAIN_TENSOR", + "weight" : 1.0, + "variable_data" : {} + }"""), + Kratos.Parameters("""{ + "type" : "strain_sensor", + "name" : "strain_xz_1", + "value" : 0, + "location" : [0.5, 0.5, 0.5], + "strain_type" : "strain_xz", + "strain_variable": "GREEN_LAGRANGE_STRAIN_TENSOR", + "weight" : 1.0, + "variable_data" : {} + }"""), + Kratos.Parameters("""{ + "type" : "strain_sensor", + "name" : "strain_xz_2", + "value" : 0, + "location" : [1.5, 0.5, 0.5], + "strain_type" : "strain_xz", + "strain_variable": "GREEN_LAGRANGE_STRAIN_TENSOR", + "weight" : 1.0, + "variable_data" : {} + }"""), + Kratos.Parameters("""{ + "type" : "strain_sensor", + "name" : "strain_yz_1", + "value" : 0, + "location" : [0.5, 0.5, 0.5], + "strain_type" : "strain_yz", + "strain_variable": "GREEN_LAGRANGE_STRAIN_TENSOR", + "weight" : 1.0, + "variable_data" : {} + }"""), + Kratos.Parameters("""{ + "type" : "strain_sensor", + "name" : "strain_yz_2", + "value" : 0, + "location" : [1.5, 0.5, 0.5], + "strain_type" : "strain_yz", + "strain_variable": "GREEN_LAGRANGE_STRAIN_TENSOR", + "weight" : 1.0, + "variable_data" : {} + }"""), + Kratos.Parameters("""{ + "type" : "strain_sensor", + "name" : "strain_zz_1", + "value" : 0, + "location" : [0.5, 0.5, 0.5], + "strain_type" : "strain_zz", + "strain_variable": "GREEN_LAGRANGE_STRAIN_TENSOR", + "weight" : 1.0, + "variable_data" : {} + }"""), + Kratos.Parameters("""{ + "type" : "strain_sensor", + "name" : "strain_yz_2", + "value" : 0, + "location" : [1.5, 0.5, 0.5], + "strain_type" : "strain_zz", + "strain_variable": "GREEN_LAGRANGE_STRAIN_TENSOR", + "weight" : 1.0, + "variable_data" : {} + }""") + ] + + cls.sensors = GetSensors(cls.model_part, parameters) + cls.ref_values = [0, 6.0, 2.0, 1.0, 1.0, 3.5, 2.0, 4.5, 3.0, 2.0, 4.0, 3.0] + + def tearDown(self): + for node in self.model_part.Nodes: + node.SetSolutionStepValue(Kratos.DISPLACEMENT, node.GetValue(Kratos.DISPLACEMENT)) + + def test_CalculateValue(self): + values = [sensor.CalculateValue(self.model_part) for sensor in self.sensors] + self.assertVectorAlmostEqual(values, self.ref_values, 3) + + def test_CalculateGradient(self): + residual_matrix = Kratos.Matrix(24, 24) + response_sensitivities = Kratos.Vector() + for i, sensor in enumerate(self.sensors): + ref_value = self.ref_values[i] + delta = 1e-5 + + adjoint_element: Kratos.Element = self.adjoint_model_part.GetElement(sensor.GetValue(KratosSI.SENSOR_ELEMENT_ID)) + sensor.CalculateGradient(adjoint_element, residual_matrix, response_sensitivities, self.model_part.ProcessInfo) + for j, node in enumerate(adjoint_element.GetGeometry()): + node.SetSolutionStepValue(Kratos.DISPLACEMENT_X, node.GetSolutionStepValue(Kratos.DISPLACEMENT_X) + delta) + perturbed_value = sensor.CalculateValue(self.model_part) + sensitivity = (perturbed_value - ref_value) / delta + self.assertAlmostEqual(sensitivity, response_sensitivities[j * 3], 4) + node.SetSolutionStepValue(Kratos.DISPLACEMENT_X, node.GetSolutionStepValue(Kratos.DISPLACEMENT_X) - delta) + + node.SetSolutionStepValue(Kratos.DISPLACEMENT_Y, node.GetSolutionStepValue(Kratos.DISPLACEMENT_Y) + delta) + perturbed_value = sensor.CalculateValue(self.model_part) + sensitivity = (perturbed_value - ref_value) / delta + self.assertAlmostEqual(sensitivity, response_sensitivities[j * 3 + 1], 4) + node.SetSolutionStepValue(Kratos.DISPLACEMENT_Y, node.GetSolutionStepValue(Kratos.DISPLACEMENT_Y) - delta) + + node.SetSolutionStepValue(Kratos.DISPLACEMENT_Z, node.GetSolutionStepValue(Kratos.DISPLACEMENT_Z) + delta) + perturbed_value = sensor.CalculateValue(self.model_part) + sensitivity = (perturbed_value - ref_value) / delta + self.assertAlmostEqual(sensitivity, response_sensitivities[j * 3 + 2], 4) + node.SetSolutionStepValue(Kratos.DISPLACEMENT_Z, node.GetSolutionStepValue(Kratos.DISPLACEMENT_Z) - delta) + if __name__ == '__main__': UnitTest.main() \ No newline at end of file