strain_energy_response_function_utility.h

Go to the documentation of this file.

// KRATOS  ___|  |                   |                   |
//       \___ \  __|  __| |   |  __| __| |   |  __| _` | |
//             | |   |    |   | (    |   |   | |   (   | |
//       _____/ \__|_|   \__,_|\___|\__|\__,_|_|  \__,_|_| MECHANICS
//
//  License:         BSD License
//                   license: StructuralMechanicsApplication/license.txt
//
//  Main authors:    Baumgaertner Daniel, https://github.com/dbaumgaertner
//                   Geiser Armin, https://github.com/armingeiser
//
 
#pragma once
 
// System includes
#include <iostream>
#include <string>
 
// External includes
 
// Project includes
#include "includes/define.h"
#include "includes/kratos_parameters.h"
#include "includes/model_part.h"
#include "utilities/variable_utils.h"
#include "processes/find_nodal_neighbours_process.h"
#include "finite_difference_utility.h"
 
// ==============================================================================
 
namespace Kratos
{
 
 
 
 
 
 
 
class StrainEnergyResponseFunctionUtility
{
public:
 
    typedef array_1d<double, 3> array_3d;
 
    KRATOS_CLASS_POINTER_DEFINITION(StrainEnergyResponseFunctionUtility);
 
 
    StrainEnergyResponseFunctionUtility(ModelPart& model_part, Parameters responseSettings)
    : mrModelPart(model_part)
    {
        // Set gradient mode
        std::string gradient_mode = responseSettings["gradient_mode"].GetString();
 
        if (gradient_mode.compare("semi_analytic") == 0)
        {
            double delta = responseSettings["step_size"].GetDouble();
            mDelta = delta;
        }
        else
            KRATOS_ERROR << "Specified gradient_mode '" << gradient_mode << "' not recognized. The only option is: semi_analytic" << std::endl;
 
    }
 
    virtual ~StrainEnergyResponseFunctionUtility()
    {
    }
 
 
 
    // ==============================================================================
    void Initialize()
    {
    }
 
    // --------------------------------------------------------------------------
    double CalculateValue()
    {
        KRATOS_TRY;
 
        const ProcessInfo &CurrentProcessInfo = mrModelPart.GetProcessInfo();
        double strain_energy = 0.0;
 
        // Sum all elemental strain energy values calculated as: W_e = u_e^T K_e u_e
        for (auto& elem_i : mrModelPart.Elements())
        {
            const bool element_is_active = elem_i.IsDefined(ACTIVE) ? elem_i.Is(ACTIVE) : true;
            if(element_is_active)
            {
                Matrix LHS;
                Vector RHS;
                Vector u;
 
                // Get state solution relevant for energy calculation
                const auto& rConstElemRef = elem_i;
                rConstElemRef.GetValuesVector(u,0);
 
                elem_i.CalculateLocalSystem(LHS,RHS,CurrentProcessInfo);
 
                // Compute strain energy
                strain_energy += 0.5 * inner_prod(u,prod(LHS,u));
            }
        }
 
        return strain_energy;
 
        KRATOS_CATCH("");
    }
 
    // --------------------------------------------------------------------------
    void CalculateGradient()
    {
        KRATOS_TRY;
 
        // Formula computed in general notation:
        // \frac{dF}{dx} = \frac{\partial F}{\partial x}
        //                 - \lambda^T \frac{\partial R_S}{\partial x}
 
        // Adjoint variable:
        // \lambda       = \frac{1}{2} u^T
 
        // Correspondingly in specific notation for given response function
        // \frac{dF}{dx} = \frac{1}{2} u^T \cdot \frac{\partial f_{ext}}{\partial x}
        //                 + \frac{1}{2} u^T \cdot ( \frac{\partial f_{ext}}{\partial x} - \frac{\partial K}{\partial x} )
 
        // First gradients are initialized
        VariableUtils().SetHistoricalVariableToZero(SHAPE_SENSITIVITY, mrModelPart.Nodes());
 
 
        // Gradient calculation
        // 1st step: Calculate partial derivative of response function w.r.t. node coordinates
        // 2nd step: Calculate adjoint field
        // 3rd step: Calculate partial derivative of state equation w.r.t. node coordinates and multiply with adjoint field
 
        // Semi analytic sensitivities
        CalculateResponseDerivativePartByFiniteDifferencing();
        CalculateAdjointField();
        CalculateStateDerivativePartByFiniteDifferencing();
 
        KRATOS_CATCH("");
    }
 
    // ==============================================================================
 
 
 
 
    std::string Info() const
    {
        return "StrainEnergyResponseFunctionUtility";
    }
 
    virtual void PrintInfo(std::ostream &rOStream) const
    {
        rOStream << "StrainEnergyResponseFunctionUtility";
    }
 
    virtual void PrintData(std::ostream &rOStream) const
    {
    }
 
 
 
protected:
 
 
 
 
    // ==============================================================================
    void CalculateAdjointField()
    {
        KRATOS_TRY;
 
        // Adjoint field may be directly obtained from state solution
 
        KRATOS_CATCH("");
    }
 
    // --------------------------------------------------------------------------
    void CalculateStateDerivativePartByFiniteDifferencing()
    {
        KRATOS_TRY;
 
        // Working variables
        const ProcessInfo &CurrentProcessInfo = mrModelPart.GetProcessInfo();
 
        // Computation of: \frac{1}{2} u^T \cdot ( - \frac{\partial K}{\partial x} )
        for (auto& elem_i : mrModelPart.Elements())
        {
            const bool element_is_active = elem_i.IsDefined(ACTIVE) ? elem_i.Is(ACTIVE) : true;
            if(element_is_active)
            {
                Vector u;
                Vector lambda;
                Vector RHS;
 
                // Get state solution
                const auto& rConstElemRef = elem_i;
                rConstElemRef.GetValuesVector(u,0);
 
                // Get adjoint variables (Corresponds to 1/2*u)
                lambda = 0.5*u;
 
                // Semi-analytic computation of partial derivative of state equation w.r.t. node coordinates
                elem_i.CalculateRightHandSide(RHS, CurrentProcessInfo);
                for (auto& node_i : elem_i.GetGeometry())
                {
                    array_3d gradient_contribution(3, 0.0);
                    Vector derived_RHS = Vector(0);
 
                    // x-direction
                    FiniteDifferenceUtility::CalculateRightHandSideDerivative(elem_i, RHS, SHAPE_SENSITIVITY_X, node_i, mDelta, derived_RHS, CurrentProcessInfo);
                    gradient_contribution[0] = inner_prod(lambda, derived_RHS);
 
                    // y-direction
                    FiniteDifferenceUtility::CalculateRightHandSideDerivative(elem_i, RHS, SHAPE_SENSITIVITY_Y, node_i, mDelta, derived_RHS, CurrentProcessInfo);
                    gradient_contribution[1] = inner_prod(lambda, derived_RHS);
 
                    // z-direction
                    FiniteDifferenceUtility::CalculateRightHandSideDerivative(elem_i, RHS, SHAPE_SENSITIVITY_Z, node_i, mDelta, derived_RHS, CurrentProcessInfo);
                    gradient_contribution[2] = inner_prod(lambda, derived_RHS);
 
                    // Assemble sensitivity to node
                    noalias(node_i.FastGetSolutionStepValue(SHAPE_SENSITIVITY)) += gradient_contribution;
                }
            }
        }
 
        // Computation of \frac{1}{2} u^T \cdot ( \frac{\partial f_{ext}}{\partial x} )
        for (auto& cond_i : mrModelPart.Conditions())
        {
            //detect if the condition is active or not. If the user did not make any choice the element
            //is active by default
            const bool condition_is_active = cond_i.IsDefined(ACTIVE) ? cond_i.Is(ACTIVE) : true;
            if (condition_is_active)
            {
                Vector u;
                Vector lambda;
                Vector RHS;
 
                // Get adjoint variables (Corresponds to 1/2*u)
                const auto& rConstCondRef = cond_i;
                rConstCondRef.GetValuesVector(u,0);
                lambda = 0.5*u;
 
                // Semi-analytic computation of partial derivative of force vector w.r.t. node coordinates
                cond_i.CalculateRightHandSide(RHS, CurrentProcessInfo);
                for (auto& node_i : cond_i.GetGeometry())
                {
                    array_3d gradient_contribution(3, 0.0);
                    Vector perturbed_RHS = Vector(0);
 
                    // Pertubation, gradient analysis and recovery of x
                    node_i.X0() += mDelta;
                    cond_i.CalculateRightHandSide(perturbed_RHS, CurrentProcessInfo);
                    gradient_contribution[0] = inner_prod(lambda, (perturbed_RHS - RHS) / mDelta);
                    node_i.X0() -= mDelta;
 
                    // Reset pertubed vector
                    perturbed_RHS = Vector(0);
 
                    // Pertubation, gradient analysis and recovery of y
                    node_i.Y0() += mDelta;
                    cond_i.CalculateRightHandSide(perturbed_RHS, CurrentProcessInfo);
                    gradient_contribution[1] = inner_prod(lambda, (perturbed_RHS - RHS) / mDelta);
                    node_i.Y0() -= mDelta;
 
                    // Reset pertubed vector
                    perturbed_RHS = Vector(0);
 
                    // Pertubation, gradient analysis and recovery of z
                    node_i.Z0() += mDelta;
                    cond_i.CalculateRightHandSide(perturbed_RHS, CurrentProcessInfo);
                    gradient_contribution[2] = inner_prod(lambda, (perturbed_RHS - RHS) / mDelta);
                    node_i.Z0() -= mDelta;
 
                    // Assemble shape gradient to node
                    noalias(node_i.FastGetSolutionStepValue(SHAPE_SENSITIVITY)) += gradient_contribution;
                }
            }
        }
 
        KRATOS_CATCH("");
    }
 
    // --------------------------------------------------------------------------
    void CalculateResponseDerivativePartByFiniteDifferencing()
    {
        KRATOS_TRY;
 
        // Working variables
        const ProcessInfo &CurrentProcessInfo = mrModelPart.GetProcessInfo();
 
        // Computation of \frac{1}{2} u^T \cdot ( \frac{\partial f_{ext}}{\partial x} )
        for (auto& cond_i : mrModelPart.Conditions())
        {
            //detect if the condition is active or not. If the user did not make any choice the element
            //is active by default
            const bool condition_is_active = cond_i.IsDefined(ACTIVE) ? cond_i.Is(ACTIVE) : true;
            if (condition_is_active)
            {
                Vector u;
                Vector RHS;
 
                // Get state solution
                const auto& rConstCondRef = cond_i;
                rConstCondRef.GetValuesVector(u,0);
 
                // Perform finite differencing of RHS vector
                cond_i.CalculateRightHandSide(RHS, CurrentProcessInfo);
                for (auto& node_i : cond_i.GetGeometry())
                {
                    array_3d gradient_contribution(3, 0.0);
                    Vector perturbed_RHS = Vector(0);
 
                    // Pertubation, gradient analysis and recovery of x
                    node_i.X0() += mDelta;
                    cond_i.CalculateRightHandSide(perturbed_RHS, CurrentProcessInfo);
                    gradient_contribution[0] = inner_prod(0.5*u, (perturbed_RHS - RHS) / mDelta);
                    node_i.X0() -= mDelta;
 
                    // Reset pertubed vector
                    perturbed_RHS = Vector(0);
 
                    // Pertubation, gradient analysis and recovery of y
                    node_i.Y0() += mDelta;
                    cond_i.CalculateRightHandSide(perturbed_RHS, CurrentProcessInfo);
                    gradient_contribution[1] = inner_prod(0.5*u, (perturbed_RHS - RHS) / mDelta);
                    node_i.Y0() -= mDelta;
 
                    // Reset pertubed vector
                    perturbed_RHS = Vector(0);
 
                    // Pertubation, gradient analysis and recovery of z
                    node_i.Z0() += mDelta;
                    cond_i.CalculateRightHandSide(perturbed_RHS, CurrentProcessInfo);
                    gradient_contribution[2] = inner_prod(0.5*u, (perturbed_RHS - RHS) / mDelta);
                    node_i.Z0() -= mDelta;
 
                    // Assemble shape gradient to node
                    noalias(node_i.FastGetSolutionStepValue(SHAPE_SENSITIVITY)) += gradient_contribution;
                }
            }
        }
        KRATOS_CATCH("");
    }
 
    // ==============================================================================
 
 
 
 
 
private:
 
 
    ModelPart &mrModelPart;
    double mDelta;
 
 
 
 
 
 
    //      StrainEnergyResponseFunctionUtility& operator=(StrainEnergyResponseFunctionUtility const& rOther);
 
    //      StrainEnergyResponseFunctionUtility(StrainEnergyResponseFunctionUtility const& rOther);
 
 
}; // Class StrainEnergyResponseFunctionUtility
 
 
 
 
 
} // namespace Kratos.