optimization_utilities.h

Go to the documentation of this file.

// ==============================================================================
//  KratosShapeOptimizationApplication
//
//  License:         BSD License
//                   license: ShapeOptimizationApplication/license.txt
//
//  Main authors:    Baumgaertner Daniel, https://github.com/dbaumgaertner
//
// ==============================================================================
 
#ifndef OPTIMIZATION_UTILITIES_H
#define OPTIMIZATION_UTILITIES_H
 
// ------------------------------------------------------------------------------
// System includes
// ------------------------------------------------------------------------------
#include <iostream>
#include <string>
 
// ------------------------------------------------------------------------------
// Project includes
// ------------------------------------------------------------------------------
#include "includes/define.h"
#include "includes/model_part.h"
#include "spaces/ublas_space.h"
#include "linear_solvers/linear_solver.h"
 
// ==============================================================================
 
namespace Kratos
{
 
 
 
 
 
 
 
 
class KRATOS_API(SHAPE_OPTIMIZATION_APPLICATION) OptimizationUtilities
{
public:
 
    typedef array_1d<double,3> array_3d;
    typedef UblasSpace<double, Matrix, Vector> DenseSpace;
 
    KRATOS_CLASS_POINTER_DEFINITION(OptimizationUtilities);
 
 
 
 
 
    // ==============================================================================
    // General optimization operations
    // ==============================================================================
    static void ComputeControlPointUpdate(ModelPart& rModelPart, const double StepSize, const bool Normalize);
 
    // --------------------------------------------------------------------------
    static void AddFirstVariableToSecondVariable( ModelPart& rModelPart, const Variable<array_3d> &rFirstVariable, const Variable<array_3d> &rSecondVariable );
 
    // --------------------------------------------------------------------------
    static double ComputeL2NormOfNodalVariable( ModelPart& rModelPart, const Variable<array_3d> &rVariable);
 
    // --------------------------------------------------------------------------
    static double ComputeL2NormOfNodalVariable(ModelPart& rModelPart, const Variable<double> &rVariable);
 
    // --------------------------------------------------------------------------
    static double ComputeMaxNormOfNodalVariable(ModelPart& rModelPart, const Variable<array_3d> &rVariable);
 
    // --------------------------------------------------------------------------
    static double ComputeMaxNormOfNodalVariable(ModelPart& rModelPart, const Variable<double> &rVariable);
 
    // ==============================================================================
    // For running unconstrained descent methods
    // ==============================================================================
    static void ComputeSearchDirectionSteepestDescent(ModelPart& rModelPart);
 
    // ==============================================================================
    // For running penalized projection method
    // ==============================================================================
    static void ComputeProjectedSearchDirection(ModelPart& rModelPart);
 
    // --------------------------------------------------------------------------
    static double CorrectProjectedSearchDirection(ModelPart& rModelPart, const double PrevConstraintValue, const double ConstraintValue, const double CorrectionScaling, const bool IsAdaptive );
 
    // --------------------------------------------------------------------------
    static double ComputeCorrectionFactor(ModelPart& rModelPart, const double PrevConstraintValue, const double ConstraintValue, double& CorrectionScaling, const bool IsAdaptive);
 
    static void AssembleVector( ModelPart& rModelPart,
        Vector& rVector,
        const Variable<double> &rVariable);
 
    static void AssembleVector( ModelPart& rModelPart,
        Vector& rVector,
        const Variable<array_3d> &rVariable);
 
    static void AssignVectorToVariable(ModelPart& rModelPart,
        const Vector& rVector,
        const Variable<double> &rVariable);
 
    static void AssignVectorToVariable(ModelPart& rModelPart,
        const Vector& rVector,
        const Variable<array_3d> &rVariable);
 
    static void AssembleMatrix(ModelPart& rModelPart,
        Matrix& rMatrix,
        const std::vector<Variable<array_3d>*>& rVariables);
 
 
    static void CalculateProjectedSearchDirectionAndCorrection(
        Vector& rObjectiveGradient,
        Matrix& rConstraintGradients,
        Vector& rConstraintValues,
        LinearSolver<DenseSpace, DenseSpace>& rSolver,
        Vector& rProjectedSearchDirection,
        Vector& rRestoration);
 
     // ==============================================================================
    // For running relaxed gradient projection
    // ==============================================================================
 
    static void AssembleBufferMatrix( Matrix& rMatrix,
        const std::vector<double>& rVariables)
    {
        size_t VectorSize = rVariables.size();
        rMatrix = ZeroMatrix(VectorSize, VectorSize);
        IndexPartition(VectorSize).for_each([&](const int i) {rMatrix(i,i) = rVariables[i];} );
    }
 
    static void CalculateRelaxedProjectedSearchDirectionAndCorrection(
        Vector& rObjectiveGradient,
        Matrix& rConstraintGradients,
        Matrix& rRelaxationCoefficients,
        Vector& rCorrectionCoefficients,
        LinearSolver<DenseSpace, DenseSpace>& rSolver,
        Vector& rProjectedSearchDirection,
        Vector& rCorrection
        )
    {
        // local variable naming according to https://msulaiman.org/onewebmedia/GradProj_2.pdf
        Vector& nabla_f = rObjectiveGradient;
        Matrix& N = rConstraintGradients;
        Vector& s = rProjectedSearchDirection;
        Vector& c = rCorrection;
        Matrix& omega_r = rRelaxationCoefficients;
        Vector& omega_c = rCorrectionCoefficients;
 
 
        Matrix NTN = prod(trans(N), N);
        Matrix I = IdentityMatrix(N.size2());
        Matrix NTN_inv(NTN.size1(), NTN.size2());
 
        rSolver.Solve(NTN, NTN_inv, I); // solve with identity to get the inverse
 
 
        s = - (nabla_f - prod(N, Vector(prod(omega_r, Vector(prod(NTN_inv, Vector(prod(trans(N), nabla_f))))))));
 
        c = - prod(N, omega_c);
 
    }
 
    // ==============================================================================
 
 
 
 
 
 
 
 
 
 
protected:
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
private:
 
 
 
    // ==============================================================================
    // Initialized by class constructor
    // ==============================================================================
 
 
 
 
 
 
 
 
 
 
//      OptimizationUtilities& operator=(OptimizationUtilities const& rOther);
 
//      OptimizationUtilities(OptimizationUtilities const& rOther);
 
 
 
}; // Class OptimizationUtilities
 
 
 
 
 
 
 
}  // namespace Kratos.
 
#endif // OPTIMIZATION_UTILITIES_H