Inject Nodal forces into the .feb file

[gibboncode]

Inject Nodal forces into the .feb file

[Kevin-Mattheus-Moerman]

Hi @xavi87 thanks for testing GIBBON and for posting this.

If I understand you correctly you want to apply forces to a set of nodes. I've some questions and some tips to start specifying loads using GIBBON.

1. Are these forces spatially varying? For a spatially constant load see the demo: DEMO_febio_0002_beam_force, where forces are applied to the end of a bar.
2. You say the forces are applied in the local normal direction. Is this only in the initial state or do you want the forces to stay orthogonal also to the deforming surface?

Perhaps you could consider defining these forces either as tractions or as pressures. These are easily specified as spatially varying. Tractions are not "follower forces" while the pressures are (so they would remain normal to the surface in any deformed state).
The following demos use tractions/pressures on a surface: DEMO_febio_0012_disc_pressure (a spatially constant pressure), and DEMO_febio_0013_disc_pressure_varying (a spatially varying pressure).

Let me know if this helps.

[bohi87]

Hi Kevin,

First, thank you for your reply,

I'm using Gibbon and FEBio to model the formation of folds during the brain growth.

I have 3 scripts :

• In the first one (the main), I create my tetraedral mesh using tetgen, then I call 2 functions like this :

for i = 1:n
[E,v,Ft]=define_parameters(V,V0,El,Fb,dt,i);
forceMagnitude=[Ft(1,:) Ft(2,:) Ft(3,:)];
bcPrescribeList=[1:size(V,1)];
bcPrescribeMagnitude = forceMagnitude(ones(1,numel(bcPrescribeList)),:);
[V]=runFEBio(V,El,E,v,t,modelName,bcPrescribeList,bcPrescribeMagnitude);
t=t+1;
end;

here, n = number of iterations, V = nodes ,V0 = initial nodes ,El = tetraedral elements, Fb = faces, dt = step of growth, (E,v) = Young mod & Poisson ration for the Neo-Hookean material, Ft = forces that I want to apply to my nodes. These forces are varying spatially like the materials.

So, I get the output parameters of the first function, and I inject them into the function (runFEBio) which will generate a model with spatially varying material parameters (the same in the DEMO_spatially_varying_material_parameters), but when at the Load Section, I don't know how to aplly the forces that I already calculate from the first function.

I attach a PNG picture to explain the interaction between the two fuctions.
PS : I have spoken to Steve Maas many times by emails, he has helped me a lot to advance in my research project, and he advised me to contact you because you are an expert in the Gibbon toolbox :-)

Thank you in advance,
Cheers,

Amine

[Kevin-Mattheus-Moerman]

For non-spatially varying forces see the FEBio user manual (section 3.12.1) about nodal loads. There you can read that they are specified like this:

<nodal_load bc="x" node_set="set1">
   <scale lc="1">1.0</scale>
   <value>3.14</value>
</nodal_load>

As you can see if refers to a "set1" which is a node set. In GIBBON one would build it like this:

appliedForce=3.14; %The applied force
bcPrescribeList=1:10; %The nodes onto which the force is applied

%Define a node set
febio_spec.Geometry.NodeSet{1}.ATTR.name='bcPrescribeList';
febio_spec.Geometry.NodeSet{1}.node.ATTR.id=bcPrescribeList(:);

%Define nodal load
febio_spec.Loads.nodal_load{1}.ATTR.bc='x';
febio_spec.Loads.nodal_load{1}.ATTR.node_set=febio_spec.Geometry.NodeSet{1}.ATTR.name;
febio_spec.Loads.nodal_load{1}.scale.ATTR.lc=1;
febio_spec.Loads.nodal_load{1}.scale.VAL=1;
febio_spec.Loads.nodal_load{1}.value=appliedForce;

As you can see this is a constant force for the entire node set and this is a single component ('x'), but you can repeat this for all 3 force components (see DEMO_febio_0002_beam_force). Also you could loop over lots of node sets.

[Kevin-Mattheus-Moerman]

it looks like perhaps spatially varying nodal loads have not been implemented in FEBio. That would be a question for Steve Maas.
Did you see the questions I asked in my reply:

1. Do you need spatially varying nodal forces?
2. Do you need follower forces, i.e. forces that remain normal to the surface during deformation?

If yes to 1, then ask Steve Maas if it is possible. As soon as we know that I can add a DEMO to do it easily. From the above you can also see how in general XML is coded so it is easy to add.

If yes to 2, then I think you should not use nodal forces as these do not stay orthogonal. I would instead use pressures. You can easily convert your force descriptions to pressure descriptions by just dividing the magnitude of your force vectors by the local face area.

[faceAreas]=patch_area(facesMatrix,verticesMatrix);
[faceForces]=vertexToFaceMeasure(facesMatrix,nodalForceMag);
facesPressures=faceForces./faceAreas;

Once you have these spatially varying pressures check out DEMO_febio_0013_disc_pressure_varying to see how these are defined in the loads section.

[bohi87]

For question 1: I need spatially varying nodal forces, the force is not the same at each element,

For question 2: I don't need follower forces, because the forces are calculated w.r.t to the normals of the tetraedral elements, here the code :
%Apply forces to nodes Ft = zeros(3,max(El(:))); for i = 1:size(El,1) Ft(:,El(i,1)) = Ft(:,El(i,1)) + P{i}*(N1{i} + N2{i} + N3{i})'/6.0; Ft(:,El(i,2)) = Ft(:,El(i,2)) + P{i}*(N1{i} + N3{i} + N4{i})'/6.0; Ft(:,El(i,3)) = Ft(:,El(i,3)) + P{i}*(N2{i} + N3{i} + N4{i})'/6.0; Ft(:,El(i,4)) = Ft(:,El(i,4)) + P{i}*(N1{i} + N2{i} + N4{i})'/6.0; end
N1 , N2 , N3 and N4 are normals on each nodes of the tetraedre i
P{i} the stress of each tetraedre i
El = tetraedres

Thanks you
Cheers

[x17wang]

Hi Kevin,

Thank you for your reply. I want to know how to apply the specific and different nodal force[f(x,y,z)] on each node of the mesh.
I defined it like this:
for j = 1:size(V,1)
forceMagnitude=[F(1,j) F(2,j) F(3,j)];
bcPrescribeList=[1:size(V,1)]';
bcPrescribeMagnitude(j,:)=forceMagnitude(:);
end
V contains all node coordinates.
Is this correct?

Thanks you
Cheers,

Xiaoyu

[Kevin-Mattheus-Moerman]

Hi @x17wang thanks for your comments. I've tried implementing spatially varying nodal forces. This will soon be part of a demo. It seems there might be an FEBio bug whereby spatially varying forces are currently ignored (I get a "no force acting on system" warning).
I'll let you know when the demo functions and when it is added to GIBBON.

Either way, once I figure out what is wrong it is coded using something like this (i.e. it uses a NodeData data entry in the MeshData section):

        <Geometry>
        ...
		<NodeSet name="NodalLoadList">
			<node id="2"/>
			<node id="4"/>
			<node id="6"/>
			<node id="8"/>
		</NodeSet>
        </Geometry>
	<MeshData>
		<NodeData name="nodal_loads" node_set="NodalLoadList">
			<node lid="1">0.5</node>
			<node lid="2">0.2</node>
			<node lid="3">0.1</node>
			<node lid="4">0.3</node>
		</NodeData>
	</MeshData>
	<Loads>
		<nodal_load bc="z" node_set="NodalLoadList">
			<scale lc="1">1</scale>
			<value node_data="nodal_loads"/>
		</nodal_load>
	</Loads>

Which can be coded e.g. using language like:

NodalLoadList=[2 4 6 8];  %The list of nodes to apply forces to
forceData=[0.5 0.2 0.1 0.3]; %The nodal forces (z component here)

%Node sets
febio_spec.Geometry.NodeSet{1}.ATTR.name='NodalLoadList';
febio_spec.Geometry.NodeSet{1}.node.ATTR.id=NodalLoadList(:);

%MeshData
febio_spec.MeshData.NodeData{1}.ATTR.name='nodal_loads';
febio_spec.MeshData.NodeData{1}.ATTR.node_set=febio_spec.Geometry.NodeSet{1}.ATTR.name;
febio_spec.MeshData.NodeData{1}.node.VAL=forceData;
febio_spec.MeshData.NodeData{1}.node.ATTR.lid=(1:1:numel(NodalLoadList))';

%Loads
febio_spec.Loads.nodal_load{1}.ATTR.bc='z';
febio_spec.Loads.nodal_load{1}.ATTR.node_set=febio_spec.Geometry.NodeSet{1}.ATTR.name;
febio_spec.Loads.nodal_load{1}.scale.ATTR.lc=1;
febio_spec.Loads.nodal_load{1}.scale.VAL=1;
febio_spec.Loads.nodal_load{1}.value.ATTR.node_data='nodal_loads';

[x17wang]

Hi Kevin,

Thank you for your reply and thank you to develop a demo for us. I had the same problem as you ("no force acting on system"). And I also want to know if this definition (febio_spec.MeshData.NodeData{1}.node.VAL=forceData;) need to change to febio_spec.MeshData.NodeData{1}.node.VAL=forceData'; (transfer matrix forceData) in order to match with this definition: febio_spec.MeshData.NodeData{1}.node.ATTR.lid=(1:1:numel(NodalLoadList))';?

Cheers,

Xiaoyu

[Kevin-Mattheus-Moerman]

It is probably best to use columns for both, so forceData(:) or forceData' would work.
I've alerted Steve Maas (at FEBio) of this issue. I'll post here when I hear back from him.

[x17wang]

Hi Kevin,

Ok. Thank you.

Xiaoyu

[x17wang]

Hi Kevin,

I have another question which is about exporting the files .xplt, we konw that FEBio exports results files after every step. Can we control to export result files every N steps of FEBio by using GibbonCode, with N a number defined by user?

Thank you.
Cheers,

Xiaoyu

[Kevin-Mattheus-Moerman]

Hi @x17wang I fixed this. It was not an FEBio bug (thanks @belgiansteve for your help figuring that out), but a typo on my end. I updated the above comment showing how to code this. I accidentally had:
febio_spec.Loads.nodal_load{1}.value.ATTR.nodal_data='nodal_loads';
Which should be:
febio_spec.Loads.nodal_load{1}.value.ATTR.node_data='nodal_loads';
(i.e. node_data instead of nodal_data).
So if you change that in your code it should work. The demo I'll add to the next version of GIBBON (give me a couple of days) will be called DEMO_febio_0030_block_varying_nodal_forces.
Let me know if this works.

[x17wang]

Hi Kevin,

Thank you for your help.
I have another question about initial velocity, I want to save the final velocity at the end of a dynamic progress for importing into the next dynamic progress and serving as the initial velocity, how can I do using GibbonCode?

Cheers,

Xiaoyu

[Kevin-Mattheus-Moerman]

In relation to nodal forces there is now a new demo: DEMO_febio_0030_block_varying_nodal_forces.

Thanks for that follow up question. Essentially you can remember that anything that can be done with FEBio can be done through GIBBON. First you need to figure out how the FEBio side of things works. This leads to understanding what the FEB input file (and XML file) should look like and what entries it should have. Then in GIBBON you can code any content you like. All the of the DEMO_febio_00.... files have a code section defining a variable febio_spec which forms the FEB file structure. You can study those to see how to code stuff (i.e. attributes, lists, etc, that will end up in the FEB file). (see also HELP_febioStruct2xml).

If you implement this, would you like to create a DEMO for it and contribute it back to GIBBON?

Let me know if you have questions.

[x17wang]

Hi Kevin,

Thank you for your reply.

I can't find the definition of "Initial Section" in the function febStruct2febFile.m, even if I defined the part of initialization of velocity, it can't be sent to the file .feb. Do you know how can I solve this problem?

Thank you.

Xiaoyu.

[x17wang]

Hi Kevin,

I found that febioStruct2xml.m function has a part of Initial Section, but when I used this function febioStruct2xml to write "febio_spec" in .feb, I always encountered this error :

Reading file /home/x17wang/Bureau/MATLAB/GIBBON-master/data/temp/sphere7_1.feb ...FAILED!
FATAL ERROR: Failed opening input file /home/x17wang/Bureau/MATLAB/GIBBON-master/data/temp/sphere7_1.feb

And I found that the created .feb has no permission of execution. Do you know how could I solve this problem?

Thank you.

Xiaoyu

[Kevin-Mattheus-Moerman]

Is this still an issue? sorry for the late reply

[Kevin-Mattheus-Moerman]

👋 let me know if you still need help with this