fea linear elements refactor #552

jibril-b-coulibaly · 2025-02-27T00:44:07Z

This PR refactors the Spring, Bar, and some Beam elements.

This is ongoing work, we should decide when it is a good enough point to merge it

Springs

Modifications

precompute and store rest length as class attribute to avoid computation at each step (this length does not change during simulation)
refactor force calculation to limit number of length calculations and divisions

possible issues / bugs to be discussed

: the stiffness calculation for the damping does not take the timestep into account. if F = A * v, dF / dx = A / dt ?

Bars

Modifications

precompute and store axial stiffness k=EA/L to avoid computation at each step (this value does not change during simulation)
refactor force calculation to limit number of length calculations and divisions

possible issues / bugs to be discussed

: the stiffness calculation for the damping does not take the timestep into account. if F = A * v, dF / dx = A / dt ?

Euler-Bernoulli beams

Modifications

define the exact beam reference rotation using the end points and YDir vectors in the Builder, rather than using the average Y directions of the end nodes after the fact. This is more practical, but I also believe this was a bug since the average Y direction of the nodes after Gram-Schmidt has no guarantee to be colinear with user provided Y direction of the beams when the nodes pre-exist the beam.
Refactor the update rotation function: do nothing when corotation is disabled, replace rotation matrix operations with faster quaternion operations
Do once and store expensive calculations of the rotation of nodes wrt reference configuration since these never change but were performed on it every step.
refactor GetStateBlock with lambdas to avoid code duplication
rename relative rotation quaternions for clarity.
no need to normalize before sending into Gram-Schmidt, since normalization occurs there
Other things TBD

possible issues / bugs to be discussed

the Gram-Schmidt uses a char as an index, why not an int or size_t?
That index is incremented inside a while loop and used to access a 3-vector, there is a risk to access out of bounds

if for any reason the suggested vector is colinear to the first Gram-Schmidt axis, this process essentially makes the basis random, which seems wrong.

chrono/src/chrono/core/ChVector3.h

Lines 887 to 910 in 6f4d94c

    
           template <class Real> 
        
           inline void ChVector3<Real>::GetDirectionAxesAsX(ChVector3<Real>& Vx, 
        
                                                            ChVector3<Real>& Vy, 
        
                                                            ChVector3<Real>& Vz, 
        
                                                            ChVector3<Real> y_sugg) const { 
        
               Vx = *this; 
        
               bool success = Vx.Normalize(); 
        
               if (!success) 
        
                   Vx = ChVector3<Real>(1, 0, 0); 
        
               Vz.Cross(Vx, y_sugg); 
        
               success = Vz.Normalize(); 
        
               if (!success) { 
        
                   char idx = 0; 
        
                   while (!success) { 
        
                       y_sugg[idx] += 1.0; 
        
                       Vz.Cross(Vx, y_sugg); 
        
                       success = Vz.Normalize(); 
        
                       ++idx; 
        
                   } 
        
               } 
        
               Vy.Cross(Vz, Vx); 
        
           }

…ulation to limit number of expensive length calculations

…step. Refactor force calculations

This is done for simplicity but I also think there might be a bug in SetupInitial(). When the beam is created from 2 nodes, or 1 node and 1 position, the end nodes may have different Y axes, and there is no reason why this Y axis should correspond to that of the beam element since those are determined by who created these nodes.

do nothing if corotate is disabled. use cheaper quaternion representation of rotation rather than matrix multiplication

…de duplication

tasora · 2025-02-27T14:08:57Z

Hi - question: in ChElementSpring, you changed
ChVector3d int_forceA = dir * internal_force_local;
ChVector3d int_forceB = -dir * internal_force_local;
into
ChVector3d int_forceA = dPos * internal_force_local;
ChVector3d int_forceB = -dPos * internal_force_local;
but I think that the correct form is the original one, where dir is a unit length vector representing direction, whereas dPos is not unit length.
Or have I misunderstood?

tasora · 2025-02-27T14:11:20Z

In ChElementBar, same question... you compute the two forces via
ChVector3d int_forceA = dPos * internal_force_local;
ChVector3d int_forceB = -dPos * internal_force_local;
but it is different from the original
ChVector3d int_forceA = dir * internal_force_local;
ChVector3d int_forceB = -dir * internal_force_local;
Right?
A.

jibril-b-coulibaly · 2025-02-27T16:19:12Z

Hi @tasora !

The changes in the force calculation were made to avoid computing the length and scaling multiple times. In the original code below

chrono/src/chrono/fea/ChElementBar.cpp

Lines 112 to 114 in 19e970e

    
           ChVector3d dir = (nodes[1]->GetPos() - nodes[0]->GetPos()).GetNormalized(); 
        
           double L_ref = (nodes[1]->GetX0() - nodes[0]->GetX0()).Length(); 
        
           double L = (nodes[1]->GetPos() - nodes[0]->GetPos()).Length();

GetNormalized internally computes the length, and then line 114 recomputes the same length. Inside GetNormalized the vector gets scaled by 1/length, i.e., each 3 components are multiplied by 1/length.

In the refactor below

chrono/src/chrono/fea/ChElementBar.cpp

Lines 100 to 109 in f6b06bf

    
           ChVector3d dPos = nodes[1]->GetPos() - nodes[0]->GetPos(); 
        
           double L_ref = length; 
        
           double L = dPos.Length(); 
        
           double Linv = 1.0 / L; 
        
           ChVector3d dV = nodes[1]->GetPosDt() - nodes[0]->GetPosDt(); 
        
           double internal_Kforce_local = this->kstiff * (L - L_ref); 
        
           double internal_Rforce_local = (this->rdamping * this->kstiff) * Vdot(dV, dPos) * Linv; 
        
           double internal_force_local = (internal_Kforce_local + internal_Rforce_local) * Linv; 
        
           ChVector3d int_forceA = dPos * internal_force_local; 
        
           ChVector3d int_forceB = -dPos * internal_force_local;

With the non-normalized node-to-node vector dPos, we only compute the length once on 102.
the internal_force_local in 107 is scaled by Linv, so "force" here is a bit of an abuse of language but the math ends up being the same as the original.
The difference is that instead of re-scaling every component of the vector, we rescale the magnitude, doing 1 rescaling instead of 3.

We actually do 2 rescalings, with the damping on 106, so this only saves 1 multiplication by Linv compared to the original and is insignificant in terms of performance.
However, once we decide to compute the length only once by not calling GetNormalized- which does save compute - this is the more natural way to write the rest of the code in my opinion.

I ran all CORE, COLLISION, PHYSICS, Chrono joints, SMC contact in CORE module, FEA, MODAL, VEHICLE, MULTICORE, SMC contact in MULTICORE, SYNCHRONO unit tests and they all pass.

…nce for clarity. store unchanged values that were computed each step

rserban · 2025-02-27T18:33:26Z

@jibril-b-coulibaly Please note that it is not enough to run the FEA unit tests to ensure the changes you are making are correct and do not break things. Indeed, most FEA unit tests are for the ANCF formulation. Instead, you should compare results on all appropriate FEA demos.

jibril-b-coulibaly · 2025-02-27T19:15:48Z

@rserban thank you for the information.

I additionally ran demo_FEA_beamsEuler, bin/demo_FEA_beams_constr bin/demo_FEA_beams_rotor bin/demo_FEA_beams_static bin/demo_MOD_slewing_beam that looks okay.

bin/demo_FEA_beams_extrude doesn't look like it is doing much, but that behavior is the same on main so it does not seem to be affected by these changes. Looking at the code, when the extruder is updated, a new beam element is created but no orientation (referecnce quaternion) is provided. Is that a bug? This looks like it would be restricted to only work if the beams are aligned with the "World" X,Y,Z coordinates system.

chrono/src/chrono/fea/ChBuilderBeam.cpp

Lines 744 to 790 in 19e970e

    
           void ChExtruderBeamEuler::Update() { 
        
               auto node1 = beam_nodes.back(); 
        
               ChVector3d P1 = node1->GetPos(); 
        
               double d1 = (outlet.TransformPointParentToLocal(P1)).z(); 
        
               // std::cout << " d1=" << d1 << std::endl; 
        
               if (d1 >= 0) { 
        
                   double d0 = d1 - this->h; 
        
                   ChCoordsys<> C0; 
        
                   C0.rot = outlet.rot; 
        
                   C0.pos = outlet.pos + outlet.TransformPointLocalToParent(VECT_Z * d0); 
        
                   ChCoordsys<> C0_ref; 
        
                   C0_ref.rot = node1->GetX0().GetRot(); 
        
                   C0_ref.pos = node1->GetX0().GetPos() - VECT_Z * this->h; 
        
                   auto node0 = chrono_types::make_shared<ChNodeFEAxyzrot>(ChFrame<>(C0)); 
        
                   node0->SetPosDt(outlet.TransformDirectionLocalToParent(VECT_Z * this->speed)); 
        
                   node0->SetX0(ChFrame<>(C0_ref)); 
        
                   mesh->AddNode(node0); 
        
                   beam_nodes.push_back(node0); 
        
                   actuator->Initialize(node0, ground, false, ChFrame<>(C0), ChFrame<>(C0)); 
        
                   actuator->SetSpeedFunction(chrono_types::make_shared<ChFunctionConst>(this->speed)); 
        
                   actuator->SetMotionOffset(actuator->GetMotionOffset() - this->h); 
        
                   /* 
        
                   actuator->Initialize(node0, ground, false, ChFrame<>(C0), ChFrame<>(C0)); 
        
                   actuator->SetMotionFunction(chrono_types::make_shared<ChFunctionRamp>(0,this->speed)); 
        
                   actuator->SetMotionOffset(actuator->GetMotionOffset() - this->h); 
        
                   */ 
        
                   auto element = chrono_types::make_shared<ChElementBeamEuler>(); 
        
                   mesh->AddElement(element); 
        
                   beam_elems.push_back(element); 
        
                   element->SetNodes(node0, node1); 
        
                   element->SetSection(this->beam_section); 
        
                   element->SetupInitial(mysystem); 
        
                   // add collision model to node 
        
                   if (this->contactcloud) { 
        
                       contactcloud->AddNode(node0, this->contact_radius); 
        
                   } 
        
               } 
        
               mytime = mysystem->GetChTime(); 
        
           }

Looking forward, I think the purpose of unit tests is to check that functions operate the way they are supposed to, and for a refactor like the one here, I believe we should be able to rely on them. Does that mean the unit tests coverage is insufficient to test such changes and should be expanded for Euler Bernoulli beams?

If the specific functions refactored have identical behavior as the original per unit testing, there should be no issue. There is always the possibility of functions side effects and I do agree that regression tests are also necessary to see how the changes propagate to the rest of the code.
I found the demos I ran above to mostly/only provide qualitative results. And for those that provide an output, there is no automated comparison to an analytical or "gold" solution, which makes it unpractical to get good comparisons on a large number of demos (run and pipe output for different versions, and diff the files?).
Should we move towards having quantified regression tests (or naming them as such when either current unit tests and demos are actually regression tests), in addition to unit tests and demos? Something that doesn't look off by the naked eye might still be wrong.

Thank you

…otational is off

tasora · 2025-03-07T18:06:46Z

Hi Jibril Thank you for the refactoring of FEA code. I had no much time to review all the changes because I am quite busy these days. But all looks correct so far. Just keep in mind that many finite elements in Chrono were implemented many years ago, and they really need a refresh. For example I was planning to rewrite the solid elements (tetrahedrons etc. ) both to use a total lagrangian formulation and to accommodate a new design that uncouples the pde from the element, in the sense that one could apply both elasticity and thermal problems to the same element for example. I am doing something along this line in the peridynamics/ branch. Etc etc. Thanks Alessandro

…

________________________________ From: Jibril B. Coulibaly ***@***.***> Sent: Thursday, February 27, 2025 5:19:43 PM To: projectchrono/chrono ***@***.***> Cc: Alessandro TASORA ***@***.***>; Mention ***@***.***> Subject: Re: [projectchrono/chrono] fea linear elements refactor (PR #552) Hi @tasora<https://github.com/tasora> ! The changes in the force calculation were made to avoid computing the length and scaling multiple times. In the original code below https://github.com/projectchrono/chrono/blob/19e970e97778ee0064ecf03c6d2049884001353f/src/chrono/fea/ChElementBar.cpp#L112-L114 GetNormalized internally computes the length, and then line 114 recomputes the same length. Inside GetNormalized the vector gets scaled by 1/length, i.e., each 3 components are multiplied by 1/length. In the refactor below https://github.com/projectchrono/chrono/blob/f6b06bf9a1ee52fbbbeba838a5642cda76987a10/src/chrono/fea/ChElementBar.cpp#L100-L109 With the non-normalized node-to-node vector dPos, we only compute the length once on 102. the internal_force_local in 107 is scaled by Linv, so "force" here is a bit of an abuse of language but the math ends up being the same as the original. The difference is that instead of re-scaling every component of the vector, we rescale the magnitude, doing 1 rescaling instead of 3. We actually do 2 rescalings, with the damping on 106, so this only saves 1 multiplication by Linv compared to the original and is insignificant in terms of performance. However, once we decide to compute the length only once by not calling GetNormalized- which does save compute - this is the more natural way to write the rest of the code in my opinion. I ran all FEA unit tests and they all pass. — Reply to this email directly, view it on GitHub<#552 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ABHTOBFIV74HRADT4ICVM4D2R43JPAVCNFSM6AAAAABX6RKLRGVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDMOBYGQ2TKOBTGI>. You are receiving this because you were mentioned.Message ID: ***@***.***> [jibril-b-coulibaly]jibril-b-coulibaly left a comment (projectchrono/chrono#552)<#552 (comment)> Hi @tasora<https://github.com/tasora> ! The changes in the force calculation were made to avoid computing the length and scaling multiple times. In the original code below https://github.com/projectchrono/chrono/blob/19e970e97778ee0064ecf03c6d2049884001353f/src/chrono/fea/ChElementBar.cpp#L112-L114 GetNormalized internally computes the length, and then line 114 recomputes the same length. Inside GetNormalized the vector gets scaled by 1/length, i.e., each 3 components are multiplied by 1/length. In the refactor below https://github.com/projectchrono/chrono/blob/f6b06bf9a1ee52fbbbeba838a5642cda76987a10/src/chrono/fea/ChElementBar.cpp#L100-L109 With the non-normalized node-to-node vector dPos, we only compute the length once on 102. the internal_force_local in 107 is scaled by Linv, so "force" here is a bit of an abuse of language but the math ends up being the same as the original. The difference is that instead of re-scaling every component of the vector, we rescale the magnitude, doing 1 rescaling instead of 3. We actually do 2 rescalings, with the damping on 106, so this only saves 1 multiplication by Linv compared to the original and is insignificant in terms of performance. However, once we decide to compute the length only once by not calling GetNormalized- which does save compute - this is the more natural way to write the rest of the code in my opinion. I ran all FEA unit tests and they all pass. — Reply to this email directly, view it on GitHub<#552 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ABHTOBFIV74HRADT4ICVM4D2R43JPAVCNFSM6AAAAABX6RKLRGVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDMOBYGQ2TKOBTGI>. You are receiving this because you were mentioned.Message ID: ***@***.***>

jibril-b-coulibaly added 5 commits February 23, 2025 20:15

add and precompute rest length to Spring element. Refactor force calc…

c69132d

…ulation to limit number of expensive length calculations

ChElementBar: store stiffness k=EA/l instead of calculating it every …

e28b125

…step. Refactor force calculations

refactor Euler beam UpdateRotation()

57530fc

do nothing if corotate is disabled. use cheaper quaternion representation of rotation rather than matrix multiplication

refactor ChElementBeamEuler::GetStateBlock() with lambdas to avoid co…

f6b06bf

…de duplication

jibril-b-coulibaly added 2 commits February 27, 2025 11:10

forgot to remove A matrix old update

45037bd

rename quaternions for rotation from element reference to node refere…

64f9ce6

…nce for clarity. store unchanged values that were computed each step

cosmetic clean up and cheaper local displacement calculation when cor…

5154176

…otational is off

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

fea linear elements refactor #552

fea linear elements refactor #552

jibril-b-coulibaly commented Feb 27, 2025 •

edited by dariofusai93

Loading

tasora commented Feb 27, 2025

tasora commented Feb 27, 2025

jibril-b-coulibaly commented Feb 27, 2025 •

edited

Loading

rserban commented Feb 27, 2025

jibril-b-coulibaly commented Feb 27, 2025

tasora commented Mar 7, 2025 via email

	template <class Real>
	inline void ChVector3<Real>::GetDirectionAxesAsX(ChVector3<Real>& Vx,
	ChVector3<Real>& Vy,
	ChVector3<Real>& Vz,
	ChVector3<Real> y_sugg) const {
	Vx = *this;
	bool success = Vx.Normalize();
	if (!success)
	Vx = ChVector3<Real>(1, 0, 0);

	Vz.Cross(Vx, y_sugg);
	success = Vz.Normalize();
	if (!success) {
	char idx = 0;
	while (!success) {
	y_sugg[idx] += 1.0;
	Vz.Cross(Vx, y_sugg);
	success = Vz.Normalize();
	++idx;
	}
	}

	Vy.Cross(Vz, Vx);
	}

fea linear elements refactor #552

Are you sure you want to change the base?

fea linear elements refactor #552

Conversation

jibril-b-coulibaly commented Feb 27, 2025 • edited by dariofusai93 Loading

This PR refactors the Spring, Bar, and some Beam elements.

Springs

Modifications

possible issues / bugs to be discussed

Bars

Modifications

possible issues / bugs to be discussed

Euler-Bernoulli beams

Modifications

possible issues / bugs to be discussed

tasora commented Feb 27, 2025

tasora commented Feb 27, 2025

jibril-b-coulibaly commented Feb 27, 2025 • edited Loading

rserban commented Feb 27, 2025

jibril-b-coulibaly commented Feb 27, 2025

tasora commented Mar 7, 2025 via email

jibril-b-coulibaly commented Feb 27, 2025 •

edited by dariofusai93

Loading

jibril-b-coulibaly commented Feb 27, 2025 •

edited

Loading