Shape preserving constrains, i.e. nodes sliding on body surface. More...

#include <users_modules/basic_finite_elements/src/SurfaceSlidingConstrains.hpp>

Inheritance diagram for SurfaceSlidingConstrains:

Collaboration diagram for SurfaceSlidingConstrains:

Classes
struct	DriverElementOrientation
	Class implemented by user to detect face orientation. More...

struct	MyTriangleFE

struct	OpJacobian

Public Member Functions
MyTriangleFE &	getLoopFeRhs ()

MyTriangleFE &	getLoopFeLhs ()

MoFEMErrorCode	getOptions ()

	SurfaceSlidingConstrains (MoFEM::Interface &m_field, DriverElementOrientation &orientation)

virtual	~SurfaceSlidingConstrains ()=default

MoFEMErrorCode	setOperators (int tag, const std::string lagrange_multipliers_field_name, const std::string material_field_name, const double *alpha=nullptr)

MoFEMErrorCode	setOperatorsConstrainOnly (int tag, const std::string lagrange_multipliers_field_name, const std::string material_field_name)

Public Member Functions inherited from GenericSliding
	GenericSliding ()=default

virtual	~GenericSliding ()=default

Public Attributes
MoFEM::Interface &	mField

boost::shared_ptr< MyTriangleFE >	feRhsPtr

boost::shared_ptr< MyTriangleFE >	feLhsPtr

MyTriangleFE &	feRhs

MyTriangleFE &	feLhs

DriverElementOrientation &	crackFrontOrientation

double	aLpha

Detailed Description

Shape preserving constrains, i.e. nodes sliding on body surface.

Derivation and implementation of constrains preserving body surface, i.e. body shape and volume.

The idea starts form observation that body shape can be globally characterized by constant calculated as volume over its area

\[ \frac{V}{A} = C \]

Above equation expressed in integral form is

\[ \int_\Omega \textrm{d}V = C \int_\Gamma \textrm{d}S \]

where notting that,

\[ \frac{1}{3} \int_\Omega \textrm{div}[\mathbf{X}] \textrm{d}\Omega = C \int_\Gamma \textrm{d}S \]

and applying Gauss theorem we get

\[ \int_\Gamma \mathbf{X}\cdot \frac{\mathbf{N}}{\|\mathbf{N}\|} \textrm{d}\Gamma = 3C \int_\Gamma \textrm{d}S. \]

Drooping integrals on both sides, and linearizing equation, we get

\[ \frac{\mathbf{N}}{\|\mathbf{N}\|} \cdot \delta \mathbf{X} = 3C - \frac{\mathbf{N}}{\|\mathbf{N}\|}\cdot \mathbf{X} \]

where \(\delta \mathbf{X}\) is displacement sub-inctrement. Above equation is a constrain if satisfied in body shape and volume is conserved. Final form of constrain equation is

\[ \mathcal{r} = \frac{\mathbf{N}}{\|\mathbf{N}\|}\cdot \mathbf{X} - \frac{\mathbf{N_0}}{\|\mathbf{N_0}\|}\cdot \mathbf{X}_0 = \frac{\mathbf{N}}{\|\mathbf{N}\|}\cdot (\mathbf{X}-\mathbf{X}_0) \]

In the spirit of finite element method the constrain equation is multiplied by shape functions and enforce using Lagrange multiplier method

\[ \int_\Gamma \mathbf{N}^\mathsf{T}_\lambda \left( \frac{\mathbf{N}}{\|\mathbf{N}\|}\mathbf{N}_\mathbf{X}\cdot (\overline{\mathbf{X}}-\overline{\mathbf{X}}_0) \right) \|\mathbf{N}\| \textrm{d}\Gamma = \mathbf{0}. \]

Above equation is nonlinear, applying to it Taylor expansion, we can get form which can be used with Newton interactive method

\[ \begin{split} &\int_\Gamma \mathbf{N}^\mathsf{T}_\lambda \left\{ \mathbf{N}\mathbf{N}_\mathbf{X} + \left(\mathbf{X}-\mathbf{X}_0\right) \cdot \left( \textrm{Spin}\left[\frac{\partial\mathbf{X}}{\partial\xi}\right]\cdot\mathbf{B}_\eta - \textrm{Spin}\left[\frac{\partial\mathbf{X}}{\partial\eta}\right]\cdot\mathbf{B}_\xi \right) \right\} \textrm{d}\Gamma \cdot \delta \overline{\mathbf{X}}\\ = &\int_\Gamma \mathbf{N}^\mathsf{T}_\lambda \mathbf{N}\cdot(\mathbf{X}-\mathbf{X}_0) \textrm{d}\Gamma \end{split}. \]

Equation expressing forces on shape as result of constrains, as result Lagrange multiplier method have following form

\[ \begin{split} &\int_\Gamma \mathbf{N}^\mathsf{T}_\mathbf{X} \cdot \mathbf{N} \mathbf{N}_\lambda \textrm{d}\Gamma \cdot \delta\overline{\lambda}\\ + &\int_\Gamma \lambda \mathbf{N}^\mathsf{T}_\mathbf{X} \left( \textrm{Spin}\left[\frac{\partial\mathbf{X}}{\partial\xi}\right]\cdot\mathbf{B}_\eta - \textrm{Spin}\left[\frac{\partial\mathbf{X}}{\partial\eta}\right]\cdot\mathbf{B}_\xi \right) \textrm{d}\Gamma \delta\overline{\mathbf{X}}\\ = &\int_\Gamma \lambda \mathbf{N}^\mathsf{T}_\mathbf{X} \cdot \mathbf{N} \textrm{d}\Gamma \end{split} \]

Above equations are assembled into global system of equations as following

\[ \left[ \begin{array}{cc} \mathbf{K} + \mathbf{B} & \mathbf{C}^\mathsf{T} \\ \mathbf{C} + \mathbf{A} & 0 \end{array} \right] \left\{ \begin{array}{c} \delta \overline{\mathbf{X}} \\ \delta \overline{\lambda} \end{array} \right\}= \left[ \begin{array}{c} \mathbf{f} - \mathbf{C}^\mathsf{T}\overline{\lambda} \\ \overline{\mathbf{r}} \end{array} \right] \]

where

\[ \mathbf{C}= \int_\Gamma \mathbf{N}_\lambda^\mathsf{T} \mathbf{N} \cdot \mathbf{N}_\mathbf{X} \textrm{d}\Gamma, \]

\[ \mathbf{B}= \int_\Gamma \lambda (\mathbf{X}-\mathbf{X}_0)\cdot \left( \textrm{Spin}\left[\frac{\partial\mathbf{X}}{\partial\xi}\right]\cdot\mathbf{B}_\eta - \textrm{Spin}\left[\frac{\partial\mathbf{X}}{\partial\eta}\right]\cdot\mathbf{B}_\xi \right) \textrm{d}\Gamma \]

and

\[ \mathbf{A}= \int_\Gamma \mathbf{N}^\mathsf{T}_\lambda \left(\mathbf{X}-\mathbf{X}_0\right) \cdot \left( \textrm{Spin}\left[\frac{\partial\mathbf{X}}{\partial\xi}\right]\cdot\mathbf{B}_\eta - \textrm{Spin}\left[\frac{\partial\mathbf{X}}{\partial\eta}\right]\cdot\mathbf{B}_\xi \right). \]

Examples: mesh_smoothing.cpp.

Definition at line 323 of file SurfaceSlidingConstrains.hpp.

Constructor & Destructor Documentation

◆ SurfaceSlidingConstrains()

SurfaceSlidingConstrains::SurfaceSlidingConstrains	(	MoFEM::Interface &	m_field,
		DriverElementOrientation &	orientation
	)

inline

Definition at line 405 of file SurfaceSlidingConstrains.hpp.

      : mField(m_field), feRhsPtr(new MyTriangleFE(m_field)),
        feLhsPtr(new MyTriangleFE(m_field)), feRhs(*feRhsPtr), feLhs(*feLhsPtr),
        crackFrontOrientation(orientation), aLpha(1) {
 
    ierr = getOptions();
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
  }

◆ ~SurfaceSlidingConstrains()

virtual SurfaceSlidingConstrains::~SurfaceSlidingConstrains ( )

virtualdefault

Member Function Documentation

◆ getLoopFeLhs()

MyTriangleFE & SurfaceSlidingConstrains::getLoopFeLhs ( )

inline

Definition at line 387 of file SurfaceSlidingConstrains.hpp.

387{ return feLhs; }

◆ getLoopFeRhs()

MyTriangleFE & SurfaceSlidingConstrains::getLoopFeRhs ( )

inline

Definition at line 385 of file SurfaceSlidingConstrains.hpp.

385{ return feRhs; }

◆ getOptions()

MoFEMErrorCode SurfaceSlidingConstrains::getOptions ( )

inline

Definition at line 392 of file SurfaceSlidingConstrains.hpp.

                              {
    MoFEMFunctionBegin;
    ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "",
                             "Get surface sliding constrains element scaling",
                             "none");
    CHKERRQ(ierr);
    CHKERR PetscOptionsScalar("-surface_sliding_alpha", "scaling parameter", "",
                              aLpha, &aLpha, PETSC_NULL);
    ierr = PetscOptionsEnd();
    CHKERRQ(ierr);
    MoFEMFunctionReturn(0);
  }

◆ setOperators()

MoFEMErrorCode SurfaceSlidingConstrains::setOperators	(	int	tag,
		const std::string	lagrange_multipliers_field_name,
		const std::string	material_field_name,
		const double *	alpha = `nullptr`
	)

inline

Definition at line 563 of file SurfaceSlidingConstrains.hpp.

                                                             {
    MoFEMFunctionBegin;
 
    if (alpha != nullptr) {
      aLpha = *alpha;
    }
 
    boost::shared_ptr<VectorDouble> active_variables_ptr(
        new VectorDouble(3 + 9));
    boost::shared_ptr<VectorDouble> results_ptr(new VectorDouble(3 + 9));
    boost::shared_ptr<MatrixDouble> jacobian_ptr(
        new MatrixDouble(3 + 9, 3 + 9));
 
    feRhs.getOpPtrVector().clear();
    feRhs.getOpPtrVector().push_back(new OpGetActiveDofsLambda(
        lagrange_multipliers_field_name, active_variables_ptr));
    feRhs.getOpPtrVector().push_back(new OpGetActiveDofsPositions<3>(
        material_field_name, active_variables_ptr));
    feRhs.getOpPtrVector().push_back(new OpJacobian(
        tag, lagrange_multipliers_field_name, active_variables_ptr, results_ptr,
        jacobian_ptr, crackFrontOrientation, false, aLpha));
    feRhs.getOpPtrVector().push_back(
        new OpAssembleRhs<3, 9>(lagrange_multipliers_field_name, results_ptr));
    feRhs.getOpPtrVector().push_back(
        new OpAssembleRhs<3, 9>(material_field_name, results_ptr));
 
    // Adding operators to calculate the left hand side
    feLhs.getOpPtrVector().clear();
    feLhs.getOpPtrVector().push_back(new OpGetActiveDofsLambda(
        lagrange_multipliers_field_name, active_variables_ptr));
    feLhs.getOpPtrVector().push_back(new OpGetActiveDofsPositions<3>(
        material_field_name, active_variables_ptr));
    feLhs.getOpPtrVector().push_back(new OpJacobian(
        tag, lagrange_multipliers_field_name, active_variables_ptr, results_ptr,
        jacobian_ptr, crackFrontOrientation, true, aLpha));
    feLhs.getOpPtrVector().push_back(new OpAssembleLhs<3, 9>(
        lagrange_multipliers_field_name, material_field_name, jacobian_ptr));
    feLhs.getOpPtrVector().push_back(new OpAssembleLhs<3, 9>(
        material_field_name, lagrange_multipliers_field_name, jacobian_ptr));
    feLhs.getOpPtrVector().push_back(new OpAssembleLhs<3, 9>(
        material_field_name, material_field_name, jacobian_ptr));
 
    MoFEMFunctionReturn(0);
  }

◆ setOperatorsConstrainOnly()

MoFEMErrorCode SurfaceSlidingConstrains::setOperatorsConstrainOnly	(	int	tag,
		const std::string	lagrange_multipliers_field_name,
		const std::string	material_field_name
	)

inline

Definition at line 612 of file SurfaceSlidingConstrains.hpp.

                                                                 {
    MoFEMFunctionBegin;
 
    boost::shared_ptr<VectorDouble> active_variables_ptr(
        new VectorDouble(3 + 9));
    boost::shared_ptr<VectorDouble> results_ptr(new VectorDouble(3 + 9));
    boost::shared_ptr<MatrixDouble> jacobian_ptr(
        new MatrixDouble(3 + 9, 3 + 9));
 
    // Adding operators to calculate the left hand side
    feLhs.getOpPtrVector().clear();
    feLhs.getOpPtrVector().push_back(new OpGetActiveDofsLambda(
        lagrange_multipliers_field_name, active_variables_ptr));
    feLhs.getOpPtrVector().push_back(new OpGetActiveDofsPositions<3>(
        material_field_name, active_variables_ptr));
    feLhs.getOpPtrVector().push_back(new OpJacobian(
        tag, lagrange_multipliers_field_name, active_variables_ptr, results_ptr,
        jacobian_ptr, crackFrontOrientation, true, aLpha));
    feLhs.getOpPtrVector().push_back(new OpAssembleLhs<3, 9>(
        lagrange_multipliers_field_name, material_field_name, jacobian_ptr));
 
    MoFEMFunctionReturn(0);
  }