Custom BC, fixed normal displacement

**ateshian** · 06-08-2020, 10:39 AM

Hi Fluvio,

Implementing a fixed normal displacement constraint as a plugin is not going to be trivial. While I can guide you through the process, I recommend that you first try using the "sliding-elastic" contact interface with the "tension" flag set to true (1) and the friction coefficient set to 0. This will enforce exactly what you want, but at the price of creating a rigid surface that reproduces the boundary on which you want to impose the fixed normal displacement, to represent the secondary (previously called master) surface in your contact interface.

Would that be something you are interested to try first? If you still want to create your own plugin, let me know and we can discuss how it can be done.

Best,

Gerard

**Fluvio** · 06-08-2020, 02:22 PM

Gerard,

We had discussed the use of the sliding elastic BC previously. I have yet to try your suggestion on doubling-up the sliding-elastic contact to enable friction with the tension.
We are still having some issues fine-tuning the penalty and augmented Lagrangian parameters to achieve convergence.

I came across a post in which Steve suggested creating a BC plug-in for a different application.
My thought was we could try to solve the BC and friction problems in a single plug-in.

In addition to complexity, is there another reason why you would not recommend the plug-in?
I am thinking we could try both and validate results using the sliding-elastic contact.

Thank you!

**ateshian** · 06-08-2020, 02:29 PM

Hi Fluvio,

The only issue is complexity. Since the surface normals may re-orient themselves with deformation, the fixed normal displacement BC would need to be implemented as a nonlinear surface constraint. This means that you need to first evaluate, from theory, the residual force and stiffness matrix for this constraint. Then you need to discretize the resulting surface integrals and implement them in the code. There are some existing surface constraints (and many sliding contact interfaces) in FEBio that could serve as a template for you, so you would need to familiarize yourself with those.

If you are interested in pursuing this, I can provide some guidance. Let me know.

Best,

Gerard

**Fluvio** · 06-08-2020, 03:06 PM

Gerard,

We are interested!
Let me know how to go about it and if we have to review some particular documentation or examples first.

Thank you,

**ateshian** · 06-08-2020, 03:54 PM

Hi Fluvio,

That's great, I am glad you are interested. I have put together a theoretical derivation (FEFixedNormalDisplacement.pdf) describing the governing equations that (I believe) will satisfy this constraint. Please take a look at these notes, re-derive them to check for potential errors, and let me know if you have questions about them.

If you want to start looking for representative examples of FESurfaceConstraint classes in FEBio2 which may serve as a template for this FEFixedNormalDisplacement class, take a look at FEVolumeConstraint and FEDiscreteContact in the FEBioMech project. Also take a look at the parent classes FESurfaceConstraint and FENLConstraint. FENLConstraint is an abstract class and your FEFixedNormalDisplacement class should implement all of the functions that appear in FENLConstraint. That's where it helps to look at FEVolumeConstraint and FEDiscreteContact.

The notes that I am sharing with you are meant to assist you with implementing the functions

Code:

void Residual(FEGlobalVector& R, const FETimeInfo& tp)

and

Code:

void StiffnessMatrix(FESolver* psolver, const FETimeInfo& tp)

.

Take a look at those and let me know if you have questions.

Best,

Gerard

Attached Files

FEFixedNormalDisplacement.pdf (148.0 KB, 14 views)

**Fluvio** · 06-08-2020, 04:17 PM

Gerard,

One of my colleagues and I will start digging tomorrow.
We will be sure to follow-up with questions.

Thank you!

**tsims** · 06-10-2020, 01:50 PM

Hi Gerard!

I�m Tres and I work on Fluvio�s team, we�ve gone ahead and started trying to put together the Fixed Normal Displacement plugin, but as I�m pretty new to the framework I was hoping I could get some clarifications on a couple of questions I had as I began piecing things together.

First, as I was looking at one of the classes you recommended referencing, FEVolumeConstraint, I was curious about what elements were being accessed, the Documentation says it is surface elements, so is that each Tet that has a facet on the surface of the mesh? All tets in the selected region for the constraint?
The code in question from the FEVolumeConstraint Class:
elements_snippet.png

Also, I was wondering about when some of these functions are called with respect to the simulation, particularly the ones you mentioned us needing to write, the

Code:

void Residual(FEGlobalVector& R, const FETimeInfo& tp)

function and the

Code:

void StiffnessMatrix(FESolver* psolver, const FETimeInfo& tp)

function, are they called at the end of each time step? The start? Or something else entirely?

Any advice you can offer, or resources you can point me to is greatly appreciated! Thank you for your time and assistance!

**ateshian** · 06-10-2020, 02:12 PM

Hi Tres,

The FEVolumeConstraint class requires a user to select a surface and enforce a constraint on the volume enclosed by that surface. (This would be the same requirement for your proposed FEFixedNormalDisplacement class.) So the information read from the model file into FEBio is a set of surface elements stored in an FESurface class object (available in the FECore project). One can retrieve each surface element of an FESurface object as an FESurfaceElement (also in the FECore project). FESurfaceElement has functions that return the nodal coordinates and the shape functions for that surface element, among other useful properties. You will need to familiarize yourself with that class because you will need those functions for your application.

Usually, in a Newton or quasi-Newton solution scheme, multiple iterations are needed to solve for the displacement at each time step. Each iteration solves for the displacement increment \Delta \mathbf{u}, which is used cumulatively to obtain the displacement \mathbf{u} for that time step.

The Residual function is called at the start of every iteration of the quasi-Newton (or full Newton) solver. It evaluates the contribution of FEVolumeConstraint to the global load vector. This is the function where one evaluates the virtual work expressions in equations (12) and (13) of the notes I shared with you.

The StiffnessMatrix function is called whenever the stiffness matrix needs to be updated (i.e., at the start of every iteration for a full-Newton solution, or when quasi-Newton updates (BFGS or Broyden) have exceeded the user-set maximum value). It evaluates the contribution of FEVolumeConstraint to the global stiffness matrix. This is the function that would implement equations (14) and (15) in the notes.

Though I have provided these explanations, you should not have to worry when these functions are called. You only need to modify those functions to work for FEFixedNormalDisplacement, using the equations in the notes.

Let me know if you need further clarifications.

Best,

Gerard

**Fluvio** · 06-11-2020, 08:33 PM

Gerard,

We are having some issues with the latest version of FEBioStudio and FEBio3.0.

FEBio3.0 is not able to read plug-ins due to having the wrong SDK installed

Code:

Default linear solver: pardiso
Failed loading plugin C:\Program Files\FEBioStudio1.0.0[BETA]\bin\plugins\SurfaceConstraint_FEBioPlugin.dll
 Reason: Invalid SDK version.

Did we build the plug-in using the wrong SDK?

Moreover, whenever we create a Volume Constraint, even one using FEBioStudio manually, we get the following error during execution;

Code:

 *************************************************************************
 *                                ERROR                                  *
 *                                                                       *
 * tag "constraint" (line 14392) : invalid value for attribute "type"    *
 *                                                                       *
 *************************************************************************

We wanted to do this so we could familiarize ourselves with the how the constraint appears on the .feb file and test our changes.

I happen to have an earlier version on FEBioStudio in my computer but it does not feature constraints... and, also, we do not seem to be able to download previous version from the FEBioStudio website...

**tsims** · 06-12-2020, 08:56 AM

Just a small update to Fluvio's post from last night, this morning I was able to load a plugin (in FEBio 2.9) and use it in place of the volume constraint, to do this the plugin had to be built using release mode, instead of debug mode in visual studio. Although we do still have issues with the default volume constraint implementation, and implementation in FEBio 3.

Thanks for your help!
~Tres

**ateshian** · 06-14-2020, 05:52 PM

Hi Tres, Fulvio,

We have not yet officially released FEBio 3.0, which is why, to the best of my knowledge, the SDKs are not yet available.

We hope to release FEBio 3.0 in July. If you would like to wait until then, that's fine. In the meantime, you can develop the plugin for FEBio 2.9?

Best,

Gerard

**Fluvio** · 06-18-2020, 09:19 AM

Gerard,

Apologies for the late reply.
Yes, we were able to set-up everything using FEBio 2.9.
We were able to test the baseline volume constraint on our model and now we move to integrating your derivations.

Here are the changes we are making to the code, please advise;

First, we are assuming that the following variables are defined in the code as;

Code:

Derivation    |     Code
-------------------------
n                   v
g1, g2              dxs, dxr
Na, Nb              both H1,2 and N1,2? depending on the function?

Within the Residual function, we are integrating Eq. 13 by;

Code:

//evaluate normal vector

// vec3d v = (dxr ^ dxs)*w[n] * p; //----> original
vec3d v = (dxr ^ dxs)*w[n] * p * m_eps; //----> new, with penalty factor

vec3d v_cross = (v ^ v); //----> new, we create a cross of the normal vector... the penalty could also go here?

// evalueate element forces
double* H = el.H(n);
for (int j = 0; j < neln; ++j)
{
	// fe[3 * j] += H[j] * v.x; //----> original
        fe[3 * j] += H[j] * v_cross.x; //----> new, evaluate only the corresponding element of the cross product vector

	fe[3 * j+1] += H[j] * v.y;
	fe[3 * j+2] += H[j] * v.z;
}

We will get to the Stiffness Matrix once we clarify Residual

Thank you!

**ateshian** · 06-18-2020, 10:47 AM

Hi Fluvio,

1. First, we are assuming that the following variables are defined in the code as;

dxr and dxs in the code correspond to g1 and g2 in the derivation, respectively. v in the code corresponds to g1 x g2 (not n) in the derivation in addition to being multiplied by the Gauss weight (for the numerical integration) and whichever additional factor was needed to enforce the volume constraint in FEVolumeConstraint. H[j] in the Residual code corresponds to Na, whereas N[i] corresponds to Na and N[j] corresponds to Nb in StiffnessMatrix.

Within the Residual function, we are integrating Eq. 13 by;

Equation (13) should be evaluated by changing the line of code in FEVolumeConstraint::Residual from

Code:

vec3d v = (dxr ^ dxs)*w[n]*p;

to

Code:

vec3d v = (dxr ^ dxs)*w[n]*eps*(u*n);

leaving the code beneath that statement unchanged.

Here, n is the unit vector along v,

Code:

vec3d n = dxr ^ dxs; n.unit();

and u is the displacement vector that needs to be evaluated beforehand (i.e., in lines of code preceding these statements).
However, you first need to extract the nodal displacement vectors un for that surface element. You can do this in the lines of code that follow the comment. "// get the nodal coordinates". So, replace

Code:

		for (int j=0; j<neln; ++j) x[j] = mesh.Node(el.m_node[j]).m_rt;

with

Code:

        for (int j=0; j<neln; ++j) {
            x[j] = mesh.Node(el.m_node[j]).m_rt;
            un[j] = x[j]  - mesh.Node(el.m_node[j]).m_r0;
        }

where un should be declared right after the declaration for x using

Code:

vec3d un[FEElement::MAX_NODES];

Then you need to evaluate u at the integration point (in the same way that dxr and dxs are evaluated at integration points. So replace

Code:

			for (int j=0; j<neln; ++j) 
			{
				dxr += x[j]*Gr[j];
				dxs += x[j]*Gs[j];
			}

with

Code:

			for (int j=0; j<neln; ++j) 
			{
				dxr += x[j]*Gr[j];
				dxs += x[j]*Gs[j];
				u += un[j]*H[j];
			}

where u is declared and initialized as done for dxr and dxs. (Remember to move the declaration for H to precede this block of statements).

I hope this is clear enough, but don't hesitate to ask.

Best,

Gerard

**Fluvio** · 06-18-2020, 09:52 PM

Gerard,

Thank you for the guidance.
With some minor additional changes, we were able to put together the Residual function, as shown below;

Code:

void FEFixedNormalDisplacement::Residual(FEGlobalVector& R, const FETimeInfo& tp)
{
	FEMesh& mesh = *m_s.GetMesh();

	vector<double> fe;
	vector<int> lm;

	double p = m_s.m_p;

	int NE = m_s.Elements();
	vec3d x[FEElement::MAX_NODES];
	vec3d un[FEElement::MAX_NODES];													// declaration of nodal dispalcement --06/18/2020--
	for (int i = 0; i < NE; ++i)
	{
		FESurfaceElement& el = m_s.Element(i);

		// get the nodal coordinates
		int neln = el.Nodes();
		for (int j = 0; j < neln; ++j)
		{
			x[j] = mesh.Node(el.m_node[j]).m_rt;
			un[j] = x[j] - mesh.Node(el.m_node[j]).m_r0;							// extracting nodal displacements --06/18/2020--
		}

		int ndof = 3 * neln;
		fe.resize(ndof);
		zero(fe);

		double* w = el.GaussWeights();
		int nint = el.GaussPoints();
		for (int k = 0; k < nint; ++k)												// replaced iteration counter 'n' for 'k' since normal vector 'n' is defined --06/18/2020--
		{
			double* Gr = el.Gr(k);
			double* Gs = el.Gs(k);
			vec3d dxr(0, 0, 0), dxs(0, 0, 0);
			vec3d u(0, 0, 0);														// defining vector 'u' --06/18/2020--
			double* H = el.H(k);													// mode definiton of 'H' prior to evaluation of vector 'u' --06/18/2020--
			for(int j = 0; j < neln; ++j)
			{
				dxr += x[j] * Gr[j];
				dxs += x[j] * Gs[j];
				u += un[j] * H[j];													// evaluation of u at the integration point --06/18/2020--
			}

			vec3d n = (dxr ^ dxs);
			n.unit();																// definition of 'n', the unit vector along 'v' --06/18/2020--
			//evaluate normal vector
			//vec3d v = (dxr ^ dxs)*w[n] * p;										// original to volume-constraint --06/18/2020--
			vec3d v = (dxr ^ dxs)*w[k]*m_eps*(u*n);									// modified evaluation of normal vector 'v', influenced by penalty factor 'm_eps' --06/18/2020--
			
			// evalueate element forces
			//double* H = el.H(k);													// original to volume-constraint --06/18/2020--
			for (int j = 0; j < neln; ++j)
			{
				fe[3 * j] += H[j] * v.x;
				fe[3 * j+1] += H[j] * v.y;
				fe[3 * j+2] += H[j] * v.z;
			}
		}

		UnpackLM(el, lm);

		R.Assemble(el.m_node, lm, fe);
	}
}

Minor changes included;
- Changing counter 'n' for 'k' to use the vector variable 'n'
- Used 'm_eps' instead of 'eps' for the penalty

Regarding the Stiffness Matrix function, we started making similar changes, but we have some questions regarding integrating Eqn. 15;
We are still defining 'v' as;

Code:

vec3d v = (dxr ^ dxs)*w[k] * m_eps*(u*n);

And for the update of 'Kab', I am trying to group the following operators prior to evaluating the tensor;

Code:

vect3d ka = ((u*n)*I + (n^u))*(I-(n^n))*m_eps
vect3d kb = (n^n)*n*m_eps

so that;

Code:

vec3d kab;
for(int i=0; i<neln; ++i)
	for (int j = 0; j < neln; ++j)
	{
		kab = (dxr*(N[j] * Gs[i] - N[i] * Gs[j])
			- dxs*(N[j] * Gr[i] - N[i] * Gr[j]))*w[k] * 0.5*p;

		ke[3 * i][3 * j] += 0;
		ke[3 * i][3 * j + 1] += -(N[i]*ka.z*kab.z + N[i]*N[j]*kb.z) ;
		ke[3 * i][3 * j + 2] += (N[i]*ka.y*kab.y + N[i]*N[j]*kb.y);

		ke[3 * i + 1][3 * j] += (N[i]*ka.z*kab.z + N[i]*N[j]*kb.z);
		ke[3 * i + 1][3 * j + 1] += 0;
		ke[3 * i + 1][3 * j + 2] += -(N[i]*ka.x*kab.x + N[i]*N[j]*kb.x);

		ke[3 * i + 2][3 * j] += -(N[i]*ka.y*kab.y + N[i]*N[j]*kb.y);
		ke[3 * i + 2][3 * j + 1] += (N[i]*ka.x*kab.x + N[i]*N[j]*kb.x);
		ke[3 * i + 2][3 * j + 2] += 0;
	}

I am assuming

Code:

 kab = (dxr*(N[j] * Gs[i] - N[i] * Gs[j]) - dxs*(N[j] * Gr[i] - N[i] * Gr[j]))*w[k] * 0.5*p;

= A(w) tensor

Custom BC, fixed normal displacement

Custom BC, fixed normal displacement

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment