Segmentation fault in reduced order model

[Vic M]

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

[Vic M]

Debuging shows the problem in ConvertToMatrixForm()  function. Seems some of the input matrices are wrong... but which one and how to fix it ????

[Dario Mangoni]

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

[Vic M]

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);
                }
            }
        }
    }

```

[Vic M]

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

[Vic M]

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);