Hello everyone!
I'm back to my project and I'm completely lost in attempts to create
seemingly very simple system: just connect several objects together.
To find the reason I've tried to create a test program with a VERY simple
setup (see below) and surprisingly it did't work. The error looks like this
'''
Solver setup failed
LU factorization reported a problem, zero diagonal for instance
DoStepDynamics time 6065 microseconds.
Render time 1913 microseconds.
drawGrid time 143 microseconds.
'''
The code is given below.
Hope someone can explain me what I'm doing wrong and how to fix it.
Thanks in advance
________________________
Best regards,
Victor
#include "chrono/physics/ChSystemNSC.h" #include
"chrono/physics/ChLinkMate.h" #include "chrono/physics/ChLinkLock.h" #
include "chrono/physics/ChBodyEasy.h" #include
"chrono/physics/ChLinkMotorRotationSpeed.h" #include
"chrono/fea/ChElementBeamEuler.h" #include "chrono/fea/ChBuilderBeam.h" #
include "chrono/fea/ChMesh.h" #include "chrono/fea/ChElementHexaCorot_8.h"
#include "chrono_modal/ChModalAssembly.h" #include
"chrono_irrlicht/ChVisualSystemIrrlicht.h"
#include "chrono/solver/ChDirectSolverLS.h"
#include "utils/debug_harness.h"
using namespace chrono; using namespace chrono::modal; using namespace
chrono::fea; using namespace chrono::irrlicht;
int main(int argc, char* argv[]) {
double step_size = 0.05;
// turn logging on ChLogConsole consoleLogger; consoleLogger.SetCurrentLevel
(ChLog::eChLogLevel::CHMESSAGE);
SetLog(consoleLogger); // CREATE THE MODEL
// Create a Chrono::Engine physical system ChSystemNSC sys; // auto
qr_solver = chrono_types::make_shared<ChSolverSparseQR>(); auto qr_solver =
chrono_types::make_shared<ChSolverSparseLU>(); sys.SetSolver(qr_solver);
// no gravity used here sys.Set_G_acc(VNULL);
DEBUG_MESSAGE("Config started");
if(1){ auto base = chrono_types::make_shared<ChBodyEasyBox>(1, 2, 2, 200);
base->SetPos(ChVector<>(-2, 0, 0)); base->SetBodyFixed(true); sys.Add(base);
sys.ShowHierarchy(GetLog()); } DEBUG_MESSAGE("base created"); bool add_force
= false;
std::shared_ptr<ChNodeFEAxyz> inerfaceNode; std::shared_ptr<ChNodeFEAxyz>
inerfaceNodeOut;
if(1){ DEBUG_MESSAGE("payload"); // auto payload =
chrono_types::make_shared<chrono::modal::ChModalAssembly>(); //
sys.Add(payload); // payload->SetModalMode(false); auto mesh_internal =
chrono_types::make_shared<ChMesh>(); sys.Add(mesh_internal);
DEBUG_MESSAGE("adding mesh"); // payload->AddInternal(mesh_internal);
mesh_internal->SetAutomaticGravity(false);
auto mmaterial = chrono_types::make_shared<ChContinuumElastic>(); mmaterial
->Set_E(200e9); // rubber 0.01e9, steel 200e9 mmaterial->Set_v(0.4);
mmaterial->Set_RayleighDampingK(0.004); mmaterial->Set_density(1000);
mmaterial->Set_RayleighDampingM(0.0001);
DEBUG_MESSAGE("nodes");
std::shared_ptr<ChNodeFEAxyz> node1 = chrono_types::make_shared<ChNodeFEAxyz
>(); mesh_internal->AddNode(node1); node1->SetPos(ChVector<>(1,-1,-1)); std
::shared_ptr<ChNodeFEAxyz> node2 = chrono_types::make_shared<ChNodeFEAxyz
>(); mesh_internal->AddNode(node2); node2->SetPos(ChVector<>(1, 1,-1)); std
::shared_ptr<ChNodeFEAxyz> node3 = chrono_types::make_shared<ChNodeFEAxyz
>(); mesh_internal->AddNode(node3); node3->SetPos(ChVector<>(-1,1,-1)); std
::shared_ptr<ChNodeFEAxyz> node4 = chrono_types::make_shared<ChNodeFEAxyz
>(); mesh_internal->AddNode(node4); node4->SetPos(ChVector<>(-1,-1,-1));
std::shared_ptr<ChNodeFEAxyz> node5 = chrono_types::make_shared<ChNodeFEAxyz
>(); mesh_internal->AddNode(node5); node5->SetPos(ChVector<>(1,-1, 1)); std
::shared_ptr<ChNodeFEAxyz> node6 = chrono_types::make_shared<ChNodeFEAxyz
>(); mesh_internal->AddNode(node6); node6->SetPos(ChVector<>(1, 1, 1)); std
::shared_ptr<ChNodeFEAxyz> node7 = chrono_types::make_shared<ChNodeFEAxyz
>(); mesh_internal->AddNode(node7); node7->SetPos(ChVector<>(-1,1, 1)); std
::shared_ptr<ChNodeFEAxyz> node8 = chrono_types::make_shared<ChNodeFEAxyz
>(); mesh_internal->AddNode(node8); node8->SetPos(ChVector<>(-1,-1,1));
for(auto v : mesh_internal->GetNodes()) { (std::dynamic_pointer_cast<
ChNodeFEAxyz>(v))->SetPos((std::dynamic_pointer_cast<ChNodeFEAxyz>(v))->
GetPos() + ChVector<>(0,2,1)); (std::dynamic_pointer_cast<ChNodeFEAxyz>(v))
->SetMass(0.0); }
inerfaceNode = node1; inerfaceNodeOut = node8;
DEBUG_MESSAGE("element");
auto element = chrono_types::make_shared<chrono::fea::ChElementHexaCorot_8
>(); element->SetMaterial(mmaterial); element->SetNodes(node1, node2, node3,
node4, node5, node6, node7, node8);
mesh_internal->AddElement(element); DEBUG_MESSAGE("element added");
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);
if (add_force) { // Method A (simple) // The simplest way to add a force is
to use mynode->SetForce(), or to add some ChBodyLoad. // Note: this works
only for boundary nodes! inerfaceNodeOut->SetForce(ChVector<>(0, -300, 100
)); // to trigger some vibration at the free end
// Method B (advanced) // Add a force also to internal nodes that will be
removed after modal reduction. // This can be done using a callback that
will be called all times the time integrator needs it. // You will provide a
custom force writing into computed_custom_F_full vector (note: it is up to
you to use the // proper indexes) class MyCallback : public ChModalAssembly
::CustomForceFullCallback { std::shared_ptr<chrono::fea::ChNodeFEAxyz>
interfPoint; public: MyCallback(std::shared_ptr<chrono::fea::ChNodeFEAxyz> &
interfPoint_):interfPoint(interfPoint_){}; virtual void evaluate(
ChVectorDynamic<>& computed_custom_F_full, //< compute F here, size =
n_boundary_coords_w + n_internal_coords_w const ChModalAssembly& link ///<
associated modal assembly ) { // remember! assume F vector is already
properly sized, but not zeroed! computed_custom_F_full.setZero();
// just for test, assign a force to a random coordinate of F, here an
internal node // computed_custom_F_full[computed_custom_F_full.size() - 16]
= -600; ChVector<> force(0,10,30); interfPoint->SetForce(force*std::sin(link
.GetChTime()*0.1));
} };
auto my_callback = chrono_types::make_shared<MyCallback>(inerfaceNodeOut);
// payload->RegisterCallback_CustomForceFull(my_callback);
}
// ChSparseMatrix M; // payload->GetSubassemblyMassMatrix(&M); //
DEBUG_MESSAGE(M);
DEBUG_MESSAGE("payload finished"); }
sys.Set_G_acc(ChVector<>(0,0,-9.81));
sys.Setup(); sys.Update();
PRINT_VARIABLE(sys.GetNbodies()); PRINT_VARIABLE(sys.GetNbodiesFixed());
ChSparseMatrix M; sys.GetMassMatrix(&M); PRINT_VARIABLE(M);
DEBUG_MESSAGE("Config finished");
ChVisualSystemIrrlicht vis; vis.AttachSystem(&sys); vis.SetWindowSize(1024,
768); vis.SetWindowTitle("Modal reduction"); vis.Initialize(); vis.AddLogo
(); vis.ShowExplorer(true); // vis.AddSkyBox(); vis.AddCamera(ChVector<>(1,
1.3, 6), ChVector<>(3, 0, 0)); vis.AddLightWithShadow(ChVector<>(20, 20, 20
), ChVector<>(0, 0, 0), 50, 5, 50, 55); vis.AddLight(ChVector<>(-20, -20, 0
), 6, ChColor(0.6f, 1.0f, 1.0f)); vis.AddLight(ChVector<>(0, -20, -20), 6,
ChColor(0.6f, 1.0f, 1.0f));
vis.BindAll();
std::cout<<"sys.GetNdoc(): " << sys.GetNdoc() << std::endl;
std::cout<<"sys.GetNdof():
" << sys.GetNdof() << std::endl; std::cout<<"sys.GetNconstr(): " << sys.
GetNconstr() << std::endl; std::cout<<"sys.GetDOC(): " << sys.GetDOC() <<
std::endl; std::cout<<"sys.GetDOF(): " << sys.GetDOF() << std::endl; std::
cout<<"sys.GetNsysvars(Get the number of system variables (coordinates plus
the constraint multipliers, in case of quaternions).): " << sys.GetNsysvars
() << std::endl;
std::chrono::steady_clock::time_point pr_StartTime; std::chrono::
steady_clock::time_point pr_EndTime;
while (vis.Run()) { vis.BeginScene(); pr_StartTime = std::chrono::
steady_clock::now(); vis.Render(); pr_EndTime = std::chrono::steady_clock::
now(); auto Duration = std::chrono::duration_cast<std::chrono::microseconds
>(pr_EndTime-pr_StartTime); std::cout << "Render time " << Duration.count()
<< " microseconds." << std::endl;
pr_StartTime = std::chrono::steady_clock::now(); tools::drawGrid(&vis, 1, 1,
12, 12, ChCoordsys<>(ChVector<>(0, 0, 0), CH_C_PI_2, VECT_Z), ChColor(0.5f,
0.5f, 0.5f), true); pr_EndTime = std::chrono::steady_clock::now(); Duration
= std::chrono::duration_cast<std::chrono::microseconds>(pr_EndTime-
pr_StartTime); std::cout << "drawGrid time " << Duration.count() << "
microseconds." << std::endl; vis.EndScene();
pr_StartTime = std::chrono::steady_clock::now(); sys.DoStepDynamics(
step_size); pr_EndTime = std::chrono::steady_clock::now(); Duration = std::
chrono::duration_cast<std::chrono::microseconds>(pr_EndTime-pr_StartTime);
std::cout << "DoStepDynamics time " << Duration.count() << " microseconds."
<< std::endl;
ChVector<double> asmPos;
}
return 0;
}
--
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/1b96a92f-006a-4d9a-88bc-eb5855cd84b8n%40googlegroups.com.
