#include <Membrane.h>

Public Member Functions
void	applyForce (Mdouble Kn, Mdouble critDampCoeff, Mdouble Ke, Mdouble Kd, bool bendingAreaConstant)
	Calculates and applies all neccesary forces. More...

void	applyStretchForce (Mdouble Kn, Mdouble critDampCoeff)
	Apply the force due to stretching only. More...

void	applyBendForce (Mdouble Kn, Mdouble Kd, bool bendingAreaConstant)
	Apply a force due to bending. More...

void	calculateUPre (Vec3D &state, Mdouble &length, std::array< Mdouble, 2 > &faceArea)
	Calculate some prefactors for the bending penalty. More...

Mdouble	getSineHalfTheta ()
	Calculate the sine of half the bending angle. More...

void	handleParticleRemoval (unsigned int id)
	handle the partical removal More...

void	handleParticleAddition (unsigned int id, BaseParticle *p)
	Handle the particle addition. More...

void	checkActive ()
	check if the edge should calculate bending or stretch forces More...

Public Attributes
std::array< unsigned int, 2 >	vertexId

std::array< BaseParticle *, 2 >	vertex = {{nullptr}}

std::array< unsigned int, 2 >	faceVertexId

std::array< BaseParticle *, 2 >	faceVertex = {{nullptr}}

std::array< MeshTriangle *, 2 >	face = {{nullptr}}

std::array< Mdouble, 2 >	faceInitialArea

std::array< Mdouble, 8 >	uPre

Vec3D	initialState

Mdouble	initialLength

Mdouble	initialSineHalfTheta

Mdouble	effectiveMass

bool	isStretchActive

bool	isBendActive

Detailed Description

A struct containing the information needed for calculations along one edge

Member Function Documentation

◆ applyBendForce()

void Membrane::Edge::applyBendForce	(	Mdouble	Ke,
		Mdouble	Kd,
		bool	bendingAreaConstant
	)

Apply a force due to bending.

Parameters

[in]	Ke	The spring constant used for bending.
[in]	Kn	The dissipation constant used for bending.

This function applies the forces due to bending of the triangles connected by the edge. The calculation is based on the model described in doi:10.1109/SIBGRAPI.2014.20.

 {
     // Calculation according to doi:10.1109/SIBGRAPI.2014.20
     // Now bending
     // normal / face area
     if ( !isBendActive )
     {
         return;
     }
  
     unsigned int i;
  
     std::array<Vec3D, 4> u;
     std::array<Vec3D, 4> force;
  
     Vec3D initialState1 = vertex[1]->getPosition()-vertex[0]->getPosition();
     Mdouble initialLength1 = initialState1.getLength();
  
     // Do adapt when triangle changes
     std::array<Mdouble, 2> faceArea;
     if (!bendingAreaConstant)
     {
         faceArea[0] = face[0]->getArea();
         faceArea[1] = face[1]->getArea();
     }
     else
     {
         faceArea[0] = faceInitialArea[0];
         faceArea[1] = faceInitialArea[1];
     }
     calculateUPre(initialState1, initialLength1, faceArea);
     
     
     for (i=0; i<4; i++)
     {
         u[i] = face[0]->getFaceNormal() * uPre[i] +  face[1]->getFaceNormal() * uPre[4+i];
     }
     Mdouble sineHalfTheta = getSineHalfTheta();
     Mdouble prefactor = Ke*initialLength1*initialLength1 / (faceArea[0]+faceArea[1]) * ( sineHalfTheta - initialSineHalfTheta );
  
     for (i=0; i<4; i++)
     {
         force[i] = prefactor*u[i];
     }
     
     Mdouble dThetadT = Vec3D::dot(u[0], faceVertex[0]->getVelocity())
                      + Vec3D::dot(u[1], faceVertex[1]->getVelocity())
                      + Vec3D::dot(u[2], vertex[0]->getVelocity())
                      + Vec3D::dot(u[3], vertex[1]->getVelocity());
  
     prefactor = -Kd*initialLength1*dThetadT;
     for (i=0; i<4; i++)
     {
         force[i] += prefactor*u[i];
     }
  
     faceVertex[0]->addForce(force[0]);
     faceVertex[1]->addForce(force[1]);
     vertex[0]->addForce(force[2]);
     vertex[1]->addForce(force[3]);
  
 }

References Vec3D::dot(), Vec3D::getLength(), and constants::i.

Referenced by applyForce().

◆ applyForce()

void Membrane::Edge::applyForce	(	Mdouble	Kn,
		Mdouble	critDampCoeff,
		Mdouble	Ke,
		Mdouble	Kd,
		bool	bendingAreaConstant
	)

Calculates and applies all neccesary forces.

Parameters

[in]	Kn	The spring constant used for stretching.
[in]	critDampCoeff	The critical damping coefficient used for stretching.
[in]	Ke	The spring constant used for bending.
[in]	Kn	The dissipation constant used for bending.

This function applies the forces due to the mass spring system.

 {
     applyStretchForce(Kn, critDampCoeff);
     applyBendForce(Ke, Kd, bendingAreaConstant);
 }

References applyBendForce(), and applyStretchForce().

◆ applyStretchForce()

void Membrane::Edge::applyStretchForce	(	Mdouble	Kn,
		Mdouble	critDampCoeff
	)

Apply the force due to stretching only.

Parameters

[in]	Kn	The spring constant used for stretching.
[in]	critDampCoeff	The critical damping coefficient used for stretching.

This function applies the forces due to stretching of the edge

 {
     if ( !isStretchActive ) return;
     // First stretching
     Vec3D distanceVec = vertex[1]->getPosition() - vertex[0]->getPosition();
     Mdouble distance = distanceVec.getLength();
     Vec3D normal = distanceVec/distance;
     
     if (distance<0.1*initialLength)
     {
         logger(WARN, "Edge length critical: % of original length", distance/initialLength);
     }
     
     Vec3D normalDisplacement = normal*(initialLength-distance);
  
     Vec3D force = Vec3D(0,0,0);
  
     // Normal force
     force += normalDisplacement * Kn;
  
     // Damping
     Vec3D normalRelativeVelocity = normal* Vec3D::dot(normal, vertex[1]->getVelocity() - vertex[0]->getVelocity());
     Mdouble dissipationCoefficientN = 2*critDampCoeff * sqrt(effectiveMass * Kn);
     force += -dissipationCoefficientN*normalRelativeVelocity;
  
     // Normal force and Tangential Force
     vertex[0]->addForce(-force);
     vertex[1]->addForce(force);
 }

References Vec3D::dot(), Vec3D::getLength(), logger, and WARN.

Referenced by applyForce().

◆ calculateUPre()

void Membrane::Edge::calculateUPre	(	Vec3D &	state,
		Mdouble &	length,
		std::array< Mdouble, 2 > &	faceArea
	)

Calculate some prefactors for the bending penalty.

 {
     if ( face[0]==nullptr || face[1]==nullptr)
     {
         for (unsigned int i=0; i<8; i++)
         {
             uPre[i] = 0;
         }
         return;
     }
  
     uPre[0] = length/faceArea[0];
     uPre[1] = 0;
     uPre[2] =   Vec3D::dot(state, faceVertex[0]->getPosition()-vertex[1]->getPosition()) / length / faceArea[0];
     uPre[3] = - Vec3D::dot(state, faceVertex[0]->getPosition()-vertex[0]->getPosition()) / length / faceArea[0];
  
     uPre[4] = 0;
     uPre[5] = length/faceArea[1];
  
     uPre[6] = + Vec3D::dot(state, faceVertex[1]->getPosition()-vertex[1]->getPosition()) / length / faceArea[1];
     uPre[7] = - Vec3D::dot(state, faceVertex[1]->getPosition()-vertex[0]->getPosition()) / length / faceArea[1];
 }

References Vec3D::dot(), and constants::i.

◆ checkActive()

void Membrane::Edge::checkActive ( )

check if the edge should calculate bending or stretch forces

This function checks if the edge knows al required pointers to calculate and apply stretch and/or bending forces.

 {
     isStretchActive = !(!vertex[0] || !vertex[1]) && effectiveMass > 0;
     isBendActive = isStretchActive && !(!faceVertex[0] || !faceVertex[1] || !face[0] || !face[1]);
 }

Referenced by Membrane::read().

◆ getSineHalfTheta()

Mdouble Membrane::Edge::getSineHalfTheta ( )

Calculate the sine of half the bending angle.

This function calculates the sinus of half the angle between the two connected triangles.

 {
     Mdouble dotProduct = Vec3D::dot(face[0]->getFaceNormal(), face[1]->getFaceNormal());
     Mdouble sineHalfTheta;
     if (dotProduct>=1)
     {
         sineHalfTheta = 0;
     }
     else
     {
         sineHalfTheta = sqrt((1-dotProduct)/2);
     }
     
     Vec3D initialState1 = vertex[1]->getPosition()-vertex[0]->getPosition();
     if (Vec3D::dot(Vec3D::cross(face[0]->getFaceNormal(), face[1]->getFaceNormal()), initialState1) < 0)
     {
         return -sineHalfTheta;
     }
     return sineHalfTheta;
 }

References Vec3D::cross(), and Vec3D::dot().

◆ handleParticleAddition()

void Membrane::Edge::handleParticleAddition	(	unsigned int	id,
		BaseParticle *	p
	)

Handle the particle addition.

Parameters

[in]	id	The id of the added particle
[in]	p	Pointer to the added particle

This function handles the addition of a particle with the given id to the particle handler. If the particle is contained in the edge, the edge is potentially. reactivated. This function is potentially usable for MPI parallelization

 {
     unsigned int i;
     for (i=0; i<2; i++)
     {
         if (vertexId[i] == id)
         {
             vertex[i] = p;
             checkActive();
         }
         if (faceVertexId[i] == id)
         {
             faceVertex[i] = p;
             checkActive();
         }
     }
 }

References constants::i.

◆ handleParticleRemoval()

void Membrane::Edge::handleParticleRemoval ( unsigned int id )

handle the partical removal

Parameters

[in] id The id of the removed particle

This function handles the removal of a particle with the given id from the particle handler. If the particle was contained in the edge, the edge is disabled. This function is potentially usable for MPI parallelization

 {
     unsigned int i;
     for (i=0; i<2; i++)
     {
         if (vertexId[i] == id)
         {
             vertex[i] = nullptr;
             isStretchActive = false;
             isBendActive = false;
         }
         if (faceVertexId[i] == id)
         {
             faceVertex[i] = nullptr;
             isBendActive = false;
         }
     }
 }