Setting force

[Torubara-Oleg]

Hello everyone!
Can anyone help me please? I am new to FEM and trying to validate my linear elasticity problem implementation results, but I do not quite understand how to set the value of the force applied to the surface for a specific load axis.

VectorArrayCoefficient f(dim);
for (int i = 0; i < dim-1; i++)
{
f.Set(i, new ConstantCoefficient(0.0));
}
{
Vector pull_force(mesh->bdr_attributes.Max());
pull_force = 0.0;
pull_force(1) = -1.0e-2;
f.Set(dim-1, new PWConstCoefficient(pull_force));
}

There is these peace of code in ex2.cpp. Am I right if I say that vector "pull_force" defines the force which we apply to the surface and I need to set pull_force like here "pull_force(1) = -1.0e-2;" for each axis which are defined by VectorArrayCoeficient?
Is the force value that is transmitted to MFEM the force on the entire surface or on each element?

[BBBBBbruce]

Hello,
Based on my understanding, pull_force(1) is set to the material 1, it you check the material format, the "1" is the element.attribute. It would be applied to all the elements with attribute = 1. You can change it if you need to update the position forces applied in dynamic case.

The f is for the direction, e.g. x,y,z

[Torubara-Oleg]

Hi @BBBBBbruce,

Thank you for your answer

[v-dobrev]

Hi @Torubara-Oleg,

The vector coefficient f defines the force applied to all boundary elements. Using the VectorArrayCoefficient sub-class of VectorCoefficient allows us to define its components as separate scalar coefficients. The first loop sets the first dim-1 components to 0 (for all boundary attributes). The Vector pull_force defines the last component (y in 2D, z in 3D) for the different boundary attributes: the first step pull_force = 0.0 initializes all entries of the array with 0, and the line pull_force(1) = -1.0e-2 overwrites the value for boundary attribute 2 -- the convention used by PWConstCoefficient is that index k in the input Vector corresponds to boundary attribute k+1, see

mfem/fem/coefficient.hpp

Lines 105 to 108 in 26a1809

	/// Construct the constant coefficient using a vector of constants.
	/** @a c should be a vector defined by attributes, so for region with
	attribute @a i @a c[i-1] is the coefficient in that region */
	PWConstCoefficient(Vector &c)

If you want to apply nonzero forces in other components, restricted to some boundary attributes, you can use PWConstCoefficient for the components you need.

[Torubara-Oleg]

Hi @v-dobrev,

Thank you for your answer

[gokhalen]

Hi @v-dobrev, I'm having the same problem. And I am not able to understand the explanation you gave.

I'm with you until the end of the first loop which sets all components of f to zero. dim is 2.

mfem/examples/ex2.cpp

Lines 143 to 147 in ff7ed2e

	VectorArrayCoefficient f(dim);
	for (int i = 0; i < dim-1; i++)
	{
	f.Set(i, new ConstantCoefficient(0.0));
	}

The next lines are problematic for me.

mfem/examples/ex2.cpp

Lines 149 to 152 in ff7ed2e

	Vector pull_force(mesh->bdr_attributes.Max());
	pull_force = 0.0;
	pull_force(1) = -1.0e-2;
	f.Set(dim-1, new PWConstCoefficient(pull_force));

Line 149 creates pull_force of length 3 and line 150 initializes all of its components to zero.
Line 151 sets a value of -1.0e-2 for the second component.

What is happening in line 152? f(1) is being set to new PWConstCoefficient(pull_force) which is of size 3.

See:

mfem/fem/coefficient.hpp

Lines 85 to 86 in ff7ed2e

	PWConstCoefficient(Vector &c)
	{ constants.SetSize(c.Size()); constants=c; }

Now the Set method of VectorArrayCoefficient is implemented as

mfem/fem/coefficient.hpp

Line 356 in ff7ed2e

void Set(int i, Coefficient *c) { delete Coeff[i]; Coeff[i] = c; }

This seems to me that it is setting the second component of Coeff to the address returned by new PWConstCoefficient(pull_force), not the second component of pull_force. This is puzzling.

I'm sure I'm doing something wrong and that my inexperience in C++ is showing. Nevertheless, I'd be grateful for any help.

[gokhalen]

Hi @v-dobrev

I think I have a better understanding of what is going on.

mfem/examples/ex2.cpp

Lines 143 to 153 in ff7ed2e

	VectorArrayCoefficient f(dim);
	for (int i = 0; i < dim-1; i++)
	{
	f.Set(i, new ConstantCoefficient(0.0));
	}
	{
	Vector pull_force(mesh->bdr_attributes.Max());
	pull_force = 0.0;
	pull_force(1) = -1.0e-2;
	f.Set(dim-1, new PWConstCoefficient(pull_force));
	}

Here f has as many entries as boundary attributes. When the boundary attribute is 1 the integrator will intelligently choose f(0) which corresponds to new ConstantCoefficient(0.0). When the boundary attribute is 2 the integrator will intelligently choose f(1) which corresponds to new PWConstCoefficient(pull_force).

[v-dobrev]

Hi @gokhalen,

Basically, the first loop sets the first dim-1 components of the force (i.e. x-component in 2D, or x- and y-components in 3D) to be zero everywhere, regardless of the boundary attribute.

pull_force represents the last component of the force (y-component in 2D, or z-component in 3D) which is defined as a different constant on different boundary attributes. The code simply initializes this last component to zero for all boundary attributes (pull_force = 0.0;) and then overwrites the value just for boundary attribute 2 (pull_force(1) = -1.0e-2;) -- note that the value for attribute a is stored at index (a-1).

[gokhalen]

Thanks, that makes sense. I was confused.