A simple c++ shell problem

[yekose]

Hi, i am trying to make it work a little c++ code to solve a ShellMITC4 element but the system cant solve , here i am providing the code and the program output, so is there any posibility that you can see the error on my code and point out so i can proceed?

//here is the domain and nodes

Domain *theDomain = new Domain();
Node *node1 = new Node(1, 6, -5.0, 0.0, 0.0 );
Node *node2 = new Node(2, 6, 5.0, 0.0, 0.0 );
Node *node3 = new Node(3, 6, 5.0, 0.0, 10.0 );
Node *node4 = new Node(4, 6, -5.0, 0.0, 10.0 );
theDomain->addNode(node1);
theDomain->addNode(node2);
theDomain->addNode(node3);
theDomain->addNode(node4);

//constraints material and shell definitions

ElasticMembranePlateSection *theMaterial = new ElasticMembranePlateSection(1,3.0e3,0.25,1.175,1.27);
ShellMITC4 *shell = new ShellMITC4(1,1,2,3,4,*theMaterial);
theDomain->addElement(shell);
SP_Constraint *sp1 = new SP_Constraint(1,1,0.0,true);
SP_Constraint *sp2 = new SP_Constraint(2,1,0.0,true);
theDomain->addSP_Constraint(sp1);
theDomain->addSP_Constraint(sp2);

TimeSeries *theSeries = new LinearSeries();
LoadPattern *theLoadPattern = new LoadPattern(1);
theLoadPattern->setTimeSeries(theSeries);
theDomain->addLoadPattern(theLoadPattern);
Vector theLoadValues(6);
theLoadValues(0) = 0.0;
theLoadValues(1) = 0.0;
theLoadValues(2) = -10.0;
theLoadValues(3) = 0.0;
theLoadValues(4) = 0.0;
theLoadValues(5) = 0.0;
NodalLoad *theLoad = new NodalLoad(1, 3, theLoadValues);
theDomain->addNodalLoad(theLoad, 1);
theLoad = new NodalLoad(2, 4, theLoadValues);
theDomain->addNodalLoad(theLoad, 1);

//analysis part

AnalysisModel *theModel = new AnalysisModel();
EquiSolnAlgo *theSolnAlgo = new Linear();
StaticIntegrator *theIntegrator = new LoadControl(1.0, 1, 1.0, 1.0);
ConstraintHandler *theHandler = new PenaltyConstraintHandler(1.0e8,1.0e8);
RCM *theRCM = new RCM();
DOF_Numberer *theNumberer = new DOF_Numberer(*theRCM);
BandSPDLinSolver *theSolver = new BandSPDLinLapackSolver();
LinearSOE *theSOE = new BandSPDLinSOE(*theSolver);

StaticAnalysis theAnalysis(*theDomain,
*theHandler,
*theNumberer,
*theModel,
*theSolnAlgo,
*theSOE,
*theIntegrator);
int numSteps = 1;
theAnalysis.analyze(numSteps);
opserr << *theDomain;

exit(0);

////////////////////////////////////////////////////////////////////////////////////////////////
// here is the program output---------------------------------------------------------------------------------------------------

WARNING BandSPDLinLapackSolver::solve() - the LAPACK routines returned 19
WARNING Linear::solveCurrentStep() -the LinearSOE failed in solve()
StaticAnalysis::analyze() - the Algorithm failed at iteration: 0 with domain at
load factor 1
Current Domain Information
Current Time: 0
tCommitted Time: 0

NODE DATA: NumNodes: 4

Node: 1
Coordinates : -5 0 0
Disps: 0 0 0 0 0 0
unbalanced Load: 0 0 0 0 0 0
ID : 18 19 20 21 22 23


Node: 2
Coordinates : 5 0 0
Disps: 0 0 0 0 0 0
unbalanced Load: 0 0 0 0 0 0
ID : 12 13 14 15 16 17


Node: 3
Coordinates : 5 0 10
Disps: 0 0 0 0 0 0
unbalanced Load: 0 0 0 0 0 0
ID : 6 7 8 9 10 11


Node: 4
Coordinates : -5 0 10
Disps: 0 0 0 0 0 0
unbalanced Load: 0 0 0 0 0 0
ID : 0 1 2 3 4 5


ELEMENT DATA: NumEle: 1

MITC4 Non-Locking Four Node Shell
Element Number: 1
Node 1 : 1
Node 2 : 2
Node 3 : 3
Node 4 : 4
Material Information :
ElasticMembranePlateSection:
Young's Modulus E = 3000
Poisson's Ratio nu = 0.25
Thickness h = 1.175
Density rho = 1.27


SP_Constraints: numConstraints: 2
SP_Constraint: 0 Node: 1 DOF: 1 ref value: 0 current value: 0
SP_Constraint: 1 Node: 2 DOF: 1 ref value: 0 current value: 0

Pressure_Constraints: numConstraints: 0

MP_Constraints: numConstraints: 0

LOAD PATTERNS: numPatterns: 1

Load Pattern: 1
Scale Factor: 1
Linear Series: constant factor: 1
Nodal Loads:
Nodal Load: 3 load : 0 0 -10 0 0 0
Nodal Load: 4 load : 0 0 -10 0 0 0

Elemental Loads:

Single Point Constraints:

PARAMETERS: numParameters: 0

[fmk]

you have rigid body modes, i..e. nowhere near enough boundary conditions.

[yekose]

Thank your for your reply, so that means i am making a mistake while defining constraints,

SP_Constraint *sp1 = new SP_Constraint(1,1,0.0,true); So this is not the correct way to define a fixed support. Could you please advice me the correct usage to define a fixed support in the nodes so the boundary conditions will be correct to solve the soe?
thanks for advance

[fmk]

you need to add 6 SP constraints for every one of your nodes. and the dof start at 0 and go through 5 (C indexing). right now for node 1 you have the equivalent of a
"fix 1 0 1 0 0 0 0" command in the interpreter.

[yekose]

thank you so much it worked

[cobalt]

Congratulations! You have done it!
However, I even have no idea how to run the QuickMain Porject in the solution. When I compile the project, it shows so much errors.
Would you please tell me how to make QuickMain running?
Thanks!

[yekose]

It is probably because of your tcl settings. Opensees default looks for the tcl in the C:\Program Files. If you did not set it up in that directory you should configure in the project files the correct path which you have installed your tcl so while compiling if you dont get any tcl errors let me know for the other errors i ll try to help you out.

[ArdalanNejat]

I have the compelete main file of it if anyone wants it just let me know

[ArdalanNejat]

#include <stdlib.h>
#include <conio.h>
#include <OPS_Globals.h>
#include <StandardStream.h>
#include <ModelBuilder.h>
#include <SectionForceDeformation.h>
#include <ArrayOfTaggedObjects.h>

// includes for the domain classes
#include <Domain.h>
#include <Node.h>
#include <Truss.h>
#include <ElasticMaterial.h>
#include <SP_Constraint.h>
#include <LoadPattern.h>
#include <LinearSeries.h>
#include <NodalLoad.h>

// includes for the analysis classes
#include <StaticAnalysis.h>
#include <AnalysisModel.h>
#include <Linear.h>
#include <PenaltyConstraintHandler.h>
#include <DOF_Numberer.h>
#include <RCM.h>
#include <LoadControl.h>
#include <BandSPDLinSOE.h>
#include <BandSPDLinLapackSolver.h>
//#include <ElasticBeam2d.h>
///////////////////////////////////////////////////////////////////////////////
// standard C++ includes
#include <stdlib.h>
#include <iostream>
#include <conio.h>
#include <G3Globals.h>
#include <ConsoleErrorHandler.h>
#include <CrdTransf.h>
#include <ArrayOfTaggedObjects.h>
#include <MapOfTaggedObjects.h>
#include <MapOfTaggedObjectsIter.h>

// includes for the domain classes
#include <Domain.h>
#include <Node.h>
//#include <Truss.h>
//#include <Brick.h>
//MATERIAL
#include <ElasticIsotropicMaterial.h>
#include <FiberSectionRepr.h>
//#include <ElasticMaterial.h>
//a#include <ElasticIsotropic3D.h>
#include <ElasticForceBeamColumn2d.h>
#include <SP_Constraint.h>
#include <LoadPattern.h>
#include <LinearSeries.h>
#include <NodalLoad.h>

// includes for the analysis classes
#include <StaticAnalysis.h>
#include <AnalysisModel.h>
#include <Linear.h>
#include <PenaltyConstraintHandler.h>
#include <DOF_Numberer.h>
#include <RCM.h>
#include <LoadControl.h>
#include <BandSPDLinSOE.h>
#include <BandSPDLinLapackSolver.h>
#include <NDMaterial.h>
#include <MovableObject.h>
#include <TaggedObject.h>
#include <LinearCrdTransf2d.h>
#include <CorotCrdTransf2d.h>
#include <LinearCrdTransf3d.h>
//#include <ModElasticBeam2d.h>
#include <CorotCrdTransf3d.h>
#include <Beam2dUniformLoad.h>
#include <Channel.h>
#include <Beam2dPointLoad.h>
#include <ElasticSection2d.h>
#include <BeamIntegration.h>
#include <Concrete01.h>
#include <Steel01.h>
#include <Fiber.h>
#include <WideFlangeSectionIntegration.h>
#include <NDFiberSection2d.h>
#include <LinearCrdTransf2dInt.h>
#include <LegendreBeamIntegration.h>
#include <ForceBeamColumn2d.h>
#include <RCSectionIntegration.h>
#include <FiberSection2d.h>
#include <PDeltaCrdTransf2d.h>
#include <ElasticMembranePlateSection.h>
#include <ShellMITC4.h>
#include <PDeltaCrdTransf3d.h>
#include <ForceBeamColumn3d.h>
#include <ElasticBeam3d.h>
#include <Beam3dUniformLoad.h>
#include <CorotCrdTransf3d.h>
#include <LinearCrdTransf3d.h>
#include <ElasticIsotropicThreeDimensional.h>
#include <NDFiberSection3d.h>
#include <ElasticSection3d.h>
#include <Matrix.h>
#include <MP_Constraint.h>
#include <MP_ConstraintIter.h>
#include <RigidDiaphragm.h>
#include <ID.h>
#include <TrapezoidalBeamIntegration.h>
#include <NewtonCotesBeamIntegration.h>
#include <SparseGenColLinSolver.h>
#include <SuperLU.h>
#include <TransformationConstraintHandler.h>
#include <SparseGenColLinSOE.h>
#include <PlateFiberMaterial.h>
#include <MembranePlateFiberSection.h>
#include <ElasticPlateSection.h>
#include <PlainHandler.h>
#include <PlainNumberer.h>
#include <ModifiedNewton.h>
#include <CTestNormDispIncr.h>
#include <BandGenLinSolver.h>
#include <BandGenLinLapackSolver.h>
#include <ProfileSPDLinSolver.h>
#include <ProfileSPDLinDirectSolver.h>
#include <ProfileSPDLinSOE.h>



StandardStream sserr;
OPS_Stream *opserrPtr = &sserr;

int main(int argc, char **argv)
{
Domain *theDomain = new Domain();



Node *node1 = new Node(1, 6, 0, 0, 0);
Node *node2 = new Node(2, 6, 1, 0, 0);
Node *node3 = new Node(3, 6, 1, 1, 0);
Node *node4 = new Node(4, 6, 0, 1, 0);
theDomain->addNode(node1);
theDomain->addNode(node2);
theDomain->addNode(node3);
theDomain->addNode(node4);
SP_Constraint *theSP1 = new SP_Constraint(1, 0, 0.0, true);
SP_Constraint *theSP2 = new SP_Constraint(1, 1, 0.0, true);
SP_Constraint *theSP3 = new SP_Constraint(1, 2, 0.0, true);
SP_Constraint *theSP7 = new SP_Constraint(2, 0, 0.0, true);
SP_Constraint *theSP8 = new SP_Constraint(2, 1, 0.0, true);
SP_Constraint *theSP9 = new SP_Constraint(2, 2, 0.0, true);
SP_Constraint *theSP13 = new SP_Constraint(4, 0, 0.0, true);
SP_Constraint *theSP14 = new SP_Constraint(4, 1, 0.0, true);
SP_Constraint *theSP15 = new SP_Constraint(4, 2, 0.0, true);
theDomain->addSP_Constraint(theSP1);
theDomain->addSP_Constraint(theSP2);
theDomain->addSP_Constraint(theSP3);
theDomain->addSP_Constraint(theSP7);
theDomain->addSP_Constraint(theSP8);
theDomain->addSP_Constraint(theSP9);
theDomain->addSP_Constraint(theSP13);
theDomain->addSP_Constraint(theSP14);
theDomain->addSP_Constraint(theSP15);
NDMaterial *themat01 = new ElasticIsotropicMaterial(1, 2.e+11, 0.3);
NDMaterial *theplatefiber01 = new PlateFiberMaterial(4, *themat01);
SectionForceDeformation *thesec01 = new MembranePlateFiberSection(1, 0.02, *theplatefiber01);
SectionForceDeformation *theshellsection01 = new ElasticMembranePlateSection(1, 2.039E+10, 0.3, .25);
Element *theShellElement1 = new ShellMITC4(1, 1, 2, 3, 4, *thesec01);
theDomain->addElement(theShellElement1);



TimeSeries *theSeries = new LinearSeries();

LoadPattern *theLoadPattern = new LoadPattern(1,1);
theLoadPattern->setTimeSeries(theSeries);
theDomain->addLoadPattern(theLoadPattern);
Vector theLoadValues1(6);
theLoadValues1.Zero();
theLoadValues1(2) = -1000;
NodalLoad *theLoad1 = new NodalLoad(1, 3, theLoadValues1);
theDomain->addNodalLoad(theLoad1, 1);



ConstraintHandler *theHandler1 = new PlainHandler();
DOF_Numberer *theNumberer1 = new PlainNumberer();
BandGenLinSolver *theSolver1 = new BandGenLinLapackSolver();
double tol = 1e-8;
int numIter = 6, printIt = 0, normType = 2;
ConvergenceTest *theNewTest1 = new CTestNormDispIncr(tol, numIter, printIt, normType);
EquiSolnAlgo *theSolnAlgo1 = new ModifiedNewton(*theNewTest1,0);
int numIter1 = 1, dLambda = 1, minIncr = 1, maxIncr = 1;
StaticIntegrator *theStaticIntegrator1 = new LoadControl(dLambda, numIter1, minIncr, maxIncr);
AnalysisModel *theAnalysisModel1 = new AnalysisModel();
ProfileSPDLinSolver *theSolver;
theSolver = new ProfileSPDLinDirectSolver();
LinearSOE *theSOE1 = new ProfileSPDLinSOE(*theSolver);

StaticAnalysis theAnalysis1(*theDomain,
*theHandler1,
*theNumberer1,
*theAnalysisModel1,
*theSolnAlgo1,
*theSOE1,
*theStaticIntegrator1,
theNewTest1);
int numSteps = 1;
theAnalysis1.analyze(numSteps);
opserr << *theDomain;
opserr << getch();
exit(0);
}

[passenger]

Sorry, I'm having the same issue, I've tried to add six constraints to each node as described in the reply, but the problem is still not resolved

[getgrant]

The issue you're encountering appears to be related to the convergence of the analysis. The warning message "BandSPDLinLapackSolver::solve() - the LAPACK routines returned 19" suggests that the linear solver is having trouble solving the linear system of equations. Error code 19 typically indicates that the LAPACK solver encountered a singular matrix, which can happen when your model is not well-constrained or when the material properties or element formulations are causing numerical instability.

Here are some potential reasons for the convergence issue: VidMate Download VidMate https://myfiosgateway.win/

Boundary Conditions: Ensure that your boundary conditions (constraints and supports) are correctly applied. You've applied constraints on the x-displacement (DOF 1) of nodes 1 and 2, but you haven't constrained the rotation (DOF 5) or the y-displacement (DOF 2). Depending on your problem, you may need to apply additional constraints to ensure a statically determinate structure.

Element Formulation: The ShellMITC4 element is a specialized shell element, and its behavior can be sensitive to the choice of material properties and geometry. Double-check your material properties and ensure they are physically meaningful and compatible with your problem.

Load Magnitude: The applied load of -10 units might be too large for your model. You should check if the magnitude of the applied load is reasonable for the size and properties of your structure.

Integration and Element Size: Ensure that the integration scheme used for the element is appropriate for your problem. If the element's formulation requires a specific type of numerical integration (e.g., Gauss-Lobatto), make sure it is correctly specified.

Check for Geometric Nonlinearities: Depending on the problem, geometric nonlinearities may need to be considered. If your problem involves large deformations or strains, you might need to use a nonlinear analysis approach.

To diagnose and address the issue, you can try the following steps:

Check the boundary conditions and constraints to ensure they are correctly applied and that your structure is statically determinate.

Try reducing the load magnitude and see if the analysis converges with a smaller load. This can help identify if the load is the source of the convergence problem.

Experiment with different element formulations or material properties to see if they have an impact on the convergence behavior.

Consider using a more advanced nonlinear analysis method if your problem requires it.

Enable more detailed output or debugging options in your analysis software to get additional information about where the convergence is failing.

If you continue to encounter issues, consult the documentation or support resources for the finite element analysis software you are using for specific guidance on troubleshooting convergence problems with the ShellMITC4 element.

[noreenmeka]

Congratulations! You have done it!