Ok!
After some carefull consideration and application of general knowledge
about FEA software I've come up with a solution
In case someone else faces the same task I post part of my code below. Hope
it'll be helpfull
* // Create two FEM meshes: one for nodes that will be removed in modal
reduction, // the other for the nodes that will remain after modal
reduction. mesh_internal = chrono_types::make_shared<ChMesh>();
AddInternal(mesh_internal); // NOTE: MESH FOR INTERNAL NODES: USE
assembly->AddInternal() mesh_boundary =
chrono_types::make_shared<ChMesh>(); Add(mesh_boundary); // NOTE: MESH
FOR BOUNDARY NODES: USE assembly->Add()
mesh_internal->SetAutomaticGravity(false);
mesh_boundary->SetAutomaticGravity(false); std::shared_ptr<ChMesh>
mesh_loaded = chrono_types::make_shared<ChMesh>(); try{
ChMeshFileLoader::FromTetGenFile(mesh_loaded, ///< destination
mesh pathToNodesFile.c_str(),
///< name of the .node file
pathToEleFile.c_str(), ///< name of the .ele file
mmaterial, ///< material for the
created tetahedrons VNULL,
///< optional displacement of imported mesh
ChMatrix33<>(1) ///< optional
rotation/scaling of imported mesh ); } catch
(ChException myerr) { GetLog() << myerr.what(); throw
myerr; } mesh_internal->SetAutomaticGravity(false);
mesh_boundary->SetAutomaticGravity(false); //------------TODO:
internal and boundary nodes selection---------------------
std::cout<<"mesh_internal->GetNnodes() ==
"<<mesh_internal->GetNnodes()<<std::endl;
interfaceNodesIds.clear(); aggregatedNodesIds.clear();
interfaceNodesToAdd.clear(); interfaceNames.clear();
interfacePoints.clear();
if(!parse_interface_points_file(pathToInterfacePointsFile,
interfaceNames,
interfaceNodesIds,
aggregatedNodesIds)) { std::cout<<"Unable to parse
interface point file '"<<pathToInterfacePointsFile<<"'"<<std::endl;
std::string msg = "Unable to parse interface point file
'"+pathToInterfacePointsFile+"'"; throw new file_format(msg);
} if(!interfaceNodesIds.empty()) { //first
form internal and boundary meshes std::vector<int>
move2boundary; // total list of nodes to be made boundary
for(int i=0;i<interfaceNodesIds.size();++i) {
if(interfaceNodesIds[i]<0) // the point is being created not selected from
the mesh { ChVector<double>
newNodePos(0,0,0); // set it to the average of the aggregatedNodes
for(int j=0;j<aggregatedNodesIds[i].size();++j)
{
if(std::find(move2boundary.begin(),move2boundary.end(),aggregatedNodesIds[i][j])
== move2boundary.end()){
move2boundary.push_back(aggregatedNodesIds[i][j]);
} newNodePos +=
std::dynamic_pointer_cast<ChNodeFEAxyz>(mesh_loaded->GetNode(aggregatedNodesIds[i][j]))->GetX0();
} newNodePos /=
aggregatedNodesIds[i].size();
std::shared_ptr<ChNodeFEAxyzrot> interfNodePnt =
chrono_types::make_shared<ChNodeFEAxyzrot>();
interfNodePnt->SetMass(0);
interfNodePnt->SetPos(newNodePos);
mesh_boundary->AddNode(interfNodePnt);
interfacePoints.push_back(interfNodePnt); for(int
j=0;j<aggregatedNodesIds[i].size();++j) {
auto my_root = chrono_types::make_shared<ChLinkPointFrame>();
my_root->Initialize(std::dynamic_pointer_cast<ChNodeFEAxyz>(mesh_loaded->GetNode(aggregatedNodesIds[i][j])),
interfNodePnt);
AddLink(my_root); } }
else if(interfaceNodesIds[i]>=0){
std::shared_ptr<ChNodeFEAxyz> interfNode =
std::dynamic_pointer_cast<ChNodeFEAxyz>(mesh_loaded->GetNode(interfaceNodesIds[i]));
if(std::find(move2boundary.begin(),move2boundary.end(),interfaceNodesIds[i])
== move2boundary.end()){
move2boundary.push_back(interfaceNodesIds[i]); }
for(int j=0;j<aggregatedNodesIds[i].size();++j)
{
if(std::find(move2boundary.begin(),move2boundary.end(),aggregatedNodesIds[i][j])
== move2boundary.end()){
move2boundary.push_back(aggregatedNodesIds[i][j]);
} } std::shared_ptr<ChNodeFEAxyzrot>
interfNodePnt = chrono_types::make_shared<ChNodeFEAxyzrot>();
interfNodePnt->SetMass(0);
interfNodePnt->SetPos(interfNode->GetX0());
mesh_boundary->AddNode(interfNodePnt);
interfacePoints.push_back(interfNodePnt); //connect
interfNode to interface body (just like another aggregated node
auto my_root = chrono_types::make_shared<ChLinkPointFrame>();
my_root->Initialize(interfNode, interfNodePnt);
AddLink(my_root); for(int
j=0;j<aggregatedNodesIds[i].size();++j) {
std::shared_ptr<ChNodeFEAxyz> aggregatedNode =
(std::dynamic_pointer_cast<ChNodeFEAxyz>(mesh_loaded->GetNode(aggregatedNodesIds[i][j])));
auto linkInternal2Boundary =
chrono_types::make_shared<ChLinkPointFrame>();
linkInternal2Boundary->Initialize(aggregatedNode, interfNodePnt);
AddLink(linkInternal2Boundary); }
} } int nNodes = mesh_loaded->GetNnodes();
for(int i=0; i<nNodes; ++i) {
if(std::find(move2boundary.begin(),move2boundary.end(),i) ==
move2boundary.end()) mesh_internal->AddNode(
std::dynamic_pointer_cast<ChNodeFEAxyz>(mesh_loaded->GetNode(i)) );
else mesh_boundary->AddNode(
std::dynamic_pointer_cast<ChNodeFEAxyz>(mesh_loaded->GetNode(i)) );
} int nElements = mesh_loaded->GetNelements();
for(int i=0; i<nElements; ++i) {
mesh_internal->AddElement(mesh_loaded->GetElement(i)); }
} else { mesh_internal = mesh_loaded; // just
overrite it } //------------TODO: internal and boundary nodes
selection - end -------------- auto visualizeInternalA =
chrono_types::make_shared<ChVisualShapeFEA>(mesh_internal);
visualizeInternalA->SetFEMdataType(ChVisualShapeFEA::DataType::NODE_P);
visualizeInternalA->SetColorscaleMinMax(-600, 600);
visualizeInternalA->SetSmoothFaces(true);
visualizeInternalA->SetWireframe(true);
mesh_internal->AddVisualShapeFEA(visualizeInternalA); auto
visualizeInternalB =
chrono_types::make_shared<ChVisualShapeFEA>(mesh_internal);
visualizeInternalB->SetFEMglyphType(ChVisualShapeFEA::GlyphType::NODE_CSYS);
visualizeInternalB->SetFEMdataType(ChVisualShapeFEA::DataType::NONE);
visualizeInternalB->SetSymbolsThickness(0.2);
visualizeInternalB->SetSymbolsScale(0.1);
visualizeInternalB->SetZbufferHide(false);
mesh_internal->AddVisualShapeFEA(visualizeInternalB); auto
visualizeBoundaryB =
chrono_types::make_shared<ChVisualShapeFEA>(mesh_boundary);
visualizeBoundaryB->SetFEMglyphType(ChVisualShapeFEA::GlyphType::NODE_CSYS);
visualizeBoundaryB->SetFEMdataType(ChVisualShapeFEA::DataType::NONE);
visualizeBoundaryB->SetSymbolsThickness(0.4);
visualizeBoundaryB->SetSymbolsScale(4);
visualizeBoundaryB->SetZbufferHide(false);
mesh_boundary->AddVisualShapeFEA(visualizeBoundaryB);*
On Wednesday, December 13, 2023 at 5:58:15 PM UTC+3 Vic M wrote:
> Ok. Now I'm absolutely sure that the problem is with ChLink.... but have
> no clue on how to fix it.
> I've tried both AddInternalLink(my_root) and AddLink(my_root) with the
> same result:
> "SymShiftInvert: factorization failed with the given shift"
>
> If anyone knows where the problem lies It'd be very helpful
>
>
> On Tuesday, December 5, 2023 at 10:42:27 AM UTC+3 Vic M wrote:
>
>> Hi Dario,
>>
>> Thank you for fast response!
>> I'm sorry that I did give any details on the problem. I thought that it
>> was some stupid mistake of a newbie and a common solution is ready )
>> Since my first post I've done some experiments and got more understanding
>> of the problem.
>> First of all you are right "Panel" object has mistakes in its
>> implementation which I still do not know how to fix yet.
>>
>> I followed these three tutorials:
>>
>> https://github.com/projectchrono/chrono/blob/main/src/demos/modal/demo_MOD_assembly.cpp
>>
>> https://github.com/projectchrono/chrono/blob/main/src/demos/fea/demo_FEA_modal_assembly.cpp
>>
>> https://github.com/projectchrono/chrono/blob/main/src/demos/fea/demo_FEA_contacts_NSC.cpp
>>
>> From the last I've taken the loading of Tetgen mesh from my CAD model.
>> The other two - for assembly and modal reduction.
>>
>> Now to the core of my situation...
>> I need to put several "Panels" together to model their vibration. A
>> single panel contains up to 15k DoFs that is already too much for Spectra
>> and requires significant memory usage (that I've grasped experimentally but
>> not yet sure about it).
>> In order to put the panels together interface points (boundary points)
>> are needed. As described in the first link above the two types of mesh are
>> needed in the ChModalAssembly, so I did like this:
>>
>> mesh_internal = chrono_types::make_shared<ChMesh>(); AddInternal(
>> mesh_internal); // NOTE: MESH FOR INTERNAL NODES: USE
>> assembly->AddInternal()
>> mesh_boundary = chrono_types::make_shared<ChMesh>(); Add(mesh_boundary);
>> // NOTE: MESH FOR BOUNDARY NODES: USE assembly->Add()
>>
>> But then I loaded the positions and indices of vertices into the model
>> and found out that I can't easily connect ChNodeFEAxyz to ... anything
>> else (like bodies).
>> After removing all interface points the reduction goes smooth but my
>> panel is useless(
>> So I guess that the problem is in the constrains equations that are set
>> by ChLinkXXXXXXs. And I don't know yet how to add ChLinks correctly
>>
>> I put my code for adding interface points below, may be you could give me
>> a hint where I'm wrong?
>>
>> ```
>> if(!interfaceNodesIds.empty())
>> {
>> for(int i=0;i<interfaceNodesIds.size();++i)
>> {
>> if(interfaceNodesIds[i]<0) // the point is being created not
>> selected from the mesh
>> {
>> std::vector<std::shared_ptr<ChNodeFEAxyz>>
>> aggregatedNodes(aggregatedNodesIds[i].size());
>> ChVector<double> newNodePos(0,0,0); // set it to the
>> average of the aggregatedNodes
>>
>> for(int j=0;j<aggregatedNodesIds[i].size();++j)
>> {
>> aggregatedNodes[i] =
>> (std::dynamic_pointer_cast<ChNodeFEAxyz>(mesh_internal->GetNode(aggregatedNodesIds[i][j])));
>> newNodePos +=
>> std::dynamic_pointer_cast<ChNodeFEAxyz>(mesh_internal->GetNode(aggregatedNodesIds[i][j]))->GetX0();
>> }
>>
>> newNodePos /= aggregatedNodesIds[i].size(); //
>>
>> // std::shared_ptr<ChNodeFEAxyz> interfNode =
>> chrono_types::make_shared<ChNodeFEAxyz>(newNodePos);
>>
>> std::shared_ptr<ChBody> interfNodeBody =
>> chrono_types::make_shared<ChBody>();
>> interfNodeBody->SetMass(0.0001);
>>
>> interfNodeBody->SetInertia(ChMatrix33<double>::Identity()*0.00000001);
>> interfNodeBody->SetPos(newNodePos);
>> interfNodeBody->SetName(interfaceNames[i].c_str());
>> Add(interfNodeBody);
>>
>> //add body to the list of interface points for external
>> links attachment
>> interfacePoints.push_back(interfNodeBody);
>>
>>
>> //attach all the aggregated nodes to the interface body
>> (point)
>> for(int j=0;j<aggregatedNodes.size();++j)
>> {
>> mesh_boundary->AddNode(aggregatedNodes[i]);
>> auto my_root =
>> chrono_types::make_shared<ChLinkPointFrame>();
>> my_root->Initialize(aggregatedNodes[i],
>> interfNodeBody);
>> // AddInternalLink(my_root);
>> AddLink(my_root);
>> }
>> }
>> else if(interfaceNodesIds[i]>=0){
>> std::shared_ptr<ChNodeFEAxyz> interfNode =
>> std::dynamic_pointer_cast<ChNodeFEAxyz>(mesh_internal->GetNode(interfaceNodesIds[i]));
>> std::shared_ptr<ChBody> interfNodeBody =
>> chrono_types::make_shared<ChBody>();
>> interfNodeBody->SetMass(0.001);
>>
>> interfNodeBody->SetInertia(ChMatrix33<double>::Identity()*0.001);
>> interfNodeBody->SetPos(interfNode->GetX0());
>> interfNodeBody->SetName(interfaceNames[i].c_str());
>> Add(interfNodeBody);
>>
>> //add body to the list of interface points for external
>> links attachment
>> interfacePoints.push_back(interfNodeBody);
>>
>> std::vector<std::shared_ptr<ChNodeFEAxyz>>
>> aggregatedNodes(aggregatedNodesIds[i].size());
>>
>> for(int j=0;j<aggregatedNodesIds[i].size();++j)
>> {
>> aggregatedNodes[i] =
>> (std::dynamic_pointer_cast<ChNodeFEAxyz>(mesh_internal->GetNode(aggregatedNodesIds[i][j])));
>> }
>>
>>
>> // new interface point to be added to the system
>>
>> //attach the interface node to the interface body (point)
>> // mesh_boundary->AddNode(interfNode);
>> {
>> auto my_root =
>> chrono_types::make_shared<ChLinkPointFrame>();
>> my_root->Initialize(interfNode, interfNodeBody);
>> // AddLink(my_root);
>> AddInternalLink(my_root);
>>
>> }
>>
>> //attach all the aggregated nodes to the interface body
>> (point)
>> for(int j=0;j<aggregatedNodes.size();++j)
>> {
>> // mesh_boundary->AddNode(aggregatedNodes[i]);
>> auto my_root =
>> chrono_types::make_shared<ChLinkPointFrame>();
>> my_root->Initialize(aggregatedNodes[i],
>> interfNodeBody);
>> AddInternalLink(my_root);
>> // AddLink(my_root);
>> }
>> }
>> }
>> }
>> ```
>>
>>
>>
>>
>>
>>
>> On Monday, December 4, 2023 at 7:48:03 PM UTC+3 [email protected]
>> wrote:
>>
>>> Hi Victor,
>>> I'm not an expert on modal reduction, but in the meantime can you be
>>> more precise about the problem?
>>> For example:
>>>
>>> 1. you say "some of the input matrices are wrong". What "wrong"
>>> means to you? Are there missing coefficients? Is the size wrong? Does
>>> the
>>> stiffness matrix have different values to the one expected? Wa
>>> 2. are you sure that your "Pannel" object is correctly implemented?
>>> Does a normal dynamic simulation work?
>>> 3. is there any exception triggered and, if so, in which line? you
>>> say that you hit "some kind of problems". We are not really able to fix
>>> *"some
>>> kind"* of problems, but we may try to fix some *specific *problem.
>>> Is there any missing attachment in your post? How can we expect to give
>>> advices if we don't have neither your full code nor even a precise error
>>> output?
>>> 4. Did you split internal nodes from those at the boundaries? like
>>> this?
>>> auto mesh_internal = chrono_types::make_shared<ChMesh>();
>>> assembly->AddInternal(mesh_internal); // NOTE: MESH FOR
>>> INTERNAL NODES: USE assembly->AddInternal()
>>>
>>> auto mesh_boundary = chrono_types::make_shared<ChMesh>();
>>> assembly->Add(mesh_boundary); // NOTE: MESH FOR BOUNDARY NODES:
>>> USE assembly->Add()
>>> 5. which tutorial did you follow?
>>> 6. Did you try to give a look at this demo?
>>>
>>>
>>> https://github.com/projectchrono/chrono/blob/main/src/demos/modal/demo_MOD_assembly.cpp
>>>
>>> We are currently working in better modal reduction methods on a separate
>>> branch, but in the meanwhile it would be great if you provide better info
>>> about your problem.
>>>
>>> Dario
>>> Il giorno lunedì 4 dicembre 2023 alle 11:07:32 UTC+1 [email protected]
>>> ha scritto:
>>>
>>>> Debuging shows the problem in ConvertToMatrixForm() function. Seems
>>>> some of the input matrices are wrong... but which one and how to fix it
>>>> ????
>>>>
>>>> On Monday, December 4, 2023 at 12:55:56 PM UTC+3 Vic M wrote:
>>>>
>>>>> Hello everyone,
>>>>> I'm experiencing some kind of problem working with model reduction
>>>>> functionality in chrono.
>>>>> I'm using the ChModalAssembly for modeling and reduction (as
>>>>> repented in the tutorial)
>>>>> First I create the system
>>>>> ```
>>>>> ChSystemNSC sys; auto qr_solver = chrono_types::make_shared<
>>>>> ChSolverSparseQR>(); sys.SetSolver(qr_solver);
>>>>> ....
>>>>> ```
>>>>> Then creating my model object
>>>>> ```
>>>>> auto sp = chrono_types::make_shared<Pannel>();
>>>>> sys.Add(sp);
>>>>> ```
>>>>>
>>>>> The mesh is loaded from external TetGen file
>>>>> ```
>>>>> ChMeshFileLoader::FromTetGenFile(mesh_internal, ///< destination mesh
>>>>> pathToNodesFile.c_str(), ///< name of the .node file pathToEleFile.
>>>>> c_str(), ///< name of the .ele file mmaterial, ///< material for the
>>>>> created tetahedrons VNULL, ///< optional displacement of imported mesh
>>>>> ChMatrix33<>(1) ///< optional rotation/scaling of imported mesh );
>>>>> ```
>>>>> Then I do the reduction
>>>>> ```
>>>>> SetModalMode(true); SwitchModalReductionON( nModes, // The number of
>>>>> modes to retain from modal reduction, or a ChModalSolveUndamped with
>>>>> more settings ChModalDampingRayleigh(0.001, 0.005) // The damping
>>>>> model - Optional parameter: default is ChModalDampingNone(). );
>>>>> ```
>>>>>
>>>>> After the reduction the program crashes with SEGV at either
>>>>> ```
>>>>> ComputeModesDamped(10);
>>>>> ```
>>>>> or (if the above is commented out) at
>>>>>
>>>>> ```
>>>>> sys.DoStepDynamics(step_size);
>>>>> ```
>>>>> The model has 1302 nodes and 3441 elements
>>>>> It's very frustrating since I don't know where my mistake may be
>>>>>
>>>>> Thanks in advance for your help
>>>>>
>>>>> Regards,
>>>>> Victor
>>>>>
>>>>>
--
You received this message because you are subscribed to the Google Groups
"ProjectChrono" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/projectchrono/4568f2da-ba5f-44e3-8da8-8c3434feb0f7n%40googlegroups.com.
