Transient analysis help...

For developers writing C++, Fortran, Java, code who have questions or comments to make.

Moderators: silvia, selimgunay, Moderators

Post Reply
senseiwa
Posts: 8
Joined: Fri Nov 19, 2004 3:41 am
Location: Univ. Roma Tre

Transient analysis help...

Post by senseiwa »

I decided to use bbarbricks to analyze my model. I can do a static analysis but I'm quite uncomfortable with the transient one.

Can you help me? This is what I do:

Code: Select all

theSeries = new TrigSeries (0, 200, 1, 0);
GroundMotion *theMotion = new GroundMotion(0, 0, theSeries);        theLoadPattern = new UniformExcitation(*theMotion, 2, 1);
       
int dofArray[3] = {0, 1, 2};
ID *theDofID = new ID(dofArray, 3);

int *nodeArray = {0};
ID *theNodeID = new ID(nodeArray, 1);

const char *filename = "/tmp/DROP/node-0000.dat\0";
const char *data = "disp\0";

Recorder *theRecorder = new NodeRecorder(*theDofID, *theNodeID, 1, *theDomain, filename, data);

AnalysisModel     *theModel = new AnalysisModel();
ConvergenceTest *theTest = new CTestNormDispIncr (1.0e-12, 10, 0);
EquiSolnAlgo      *theSolnAlgo = new KrylovNewton (*theTest);
TransientIntegrator  *theIntegrator = new Newmark(0.5, 0.25);
ConstraintHandler *theHandler = new PenaltyConstraintHandler(1.0e12, 1.0e12);
RCM               *theRCM = new RCM();
DOF_Numberer      *theNumberer = new DOF_Numberer(*theRCM);
UmfpackGenLinSolver   *theSolver = new UmfpackGenLinSolver();
LinearSOE         *theSOE = new UmfpackGenLinSOE(*theSolver);
DirectIntegrationAnalysis theAnalysis(*theDomain, *theHandler,                                              *theNumberer, *theModel, *theSolnAlgo, *theSOE, *theIntegrator);

theAnalysis.analyze(200, 0.02);

But when executing the code, I get this output:

UmfpackGenLinSOE::setSize - n 108 nnz 5292 lVal 105840

Can someone help me? The documentation isn't really helping...
Boris
Posts: 95
Joined: Mon Jun 14, 2004 3:57 pm
Location: UC Davis

Post by Boris »

Why to you go through C++ inteface instead of regular OpenSees language interface...

If you prefere going through C++, use a debugger to see what is going on.

That printout is just giving you some statitistics for the
Umfpack solver...

Is the analysis actually runing, failing...? You do not seem to require any output...

Try putting number 1 as last argument for CTestNormDispIncr :



ConvergenceTest *theTest = new CTestNormDispIncr (1.0e-12, 10, 0);

also see this post for more info on possible printouts:

http://opensees.berkeley.edu/phpBB2/viewtopic.php?t=105


or just look into the source.

Boris
senseiwa
Posts: 8
Joined: Fri Nov 19, 2004 3:41 am
Location: Univ. Roma Tre

Post by senseiwa »

No, the simulation isn't working.

The displacement is always zero, but it's not realistic. An example:

Code: Select all

 Node: 188
        Coordinates  : 3.45919 -10.3116 9
        commitDisps: 0 0 0
        Velocities   : 0 0 0
        commitAccels: 0 0 0
         unbalanced Load: 0 0 0
        ID : 483 484 485


 Node: 189
        Coordinates  : 8.95919 -6.31564 9
        commitDisps: 0 0 0
        Velocities   : 0 0 0
        commitAccels: 0 0 0
         unbalanced Load: 0 0 0
        ID : 252 253 254
What I'm doing is this:
- create some cirlce of bricks (BbarBrick)
- create constraints: the base of the wall have 0..1 DoF set to zero: SP_Constraint(tot, i, 0, 0);
- create the simulation as you've seen

With the printid set to 1, I get a sequence of

Code: Select all

CTestNormDispIncr::test() - iteration: 1 current Norm: 0 (max: 1e-12 residual Norm: 0)
I don't see what I'm miswriting... it seems that the ground motion and loads are properly set:

Code: Select all

theSeries = new TrigSeries (0, 200, 1, 0);
GroundMotion *theMotion = new GroundMotion(0, 0, theSeries);
theLoadPattern = new UniformExcitation(*theMotion, 2, 1);
theDomain->addLoadPattern(theLoadPattern);
I'm quite a newbie of opensees... can anyone give me some hints?
fmk
Site Admin
Posts: 5884
Joined: Fri Jun 11, 2004 2:33 pm
Location: UC Berkeley
Contact:

Post by fmk »

the message that is printed out is not an error .. just a print statement someone left in to check on memory requirements .. do you have any
mass defined?
senseiwa
Posts: 8
Joined: Fri Nov 19, 2004 3:41 am
Location: Univ. Roma Tre

Post by senseiwa »

I thought it was done automagically by the bBarBrick, defining the geometry and the material property:

Code: Select all

theMaterial = new ElasticIsotropicMaterial(0, 3000, 0.3);

theElement = new BbarBrick(tot, en[0], en[1], en[2], en[3], en[4], en[5], en[6], en[7], *theMaterial);
So, if I'm wrong, and the probability is quite high :) what's the correct procedure for defining a simple experiment? I can't really find in what my code differs from the tcl source...
fmk
Site Admin
Posts: 5884
Joined: Fri Jun 11, 2004 2:33 pm
Location: UC Berkeley
Contact:

Post by fmk »

the ElasticIsotropicMaterial can take 4 args .. rho is the 4'th .. if pass it 3, rho is taken to be 0.0 .. hence 0 mass .. have a look at the material header file.
senseiwa
Posts: 8
Joined: Fri Nov 19, 2004 3:41 am
Location: Univ. Roma Tre

Post by senseiwa »

Here's a little excerpt of the modified code:

Code: Select all

theMaterial = new ElasticIsotropicMaterial(0, 3000, 0.3, 4);

Code: Select all

Current Domain Information
        Current Time: 4
tCommitted Time: 4

NODE DATA: NumNodes: 300

 Node: 0
        Coordinates  : 14 0 0
        commitDisps: 0 0 0
        Velocities   : 0 0 0
        commitAccels: 0 0 0
         unbalanced Load: 0 0 0
        ID : 222 223 224
...
 Node: 299
        Coordinates  : 8.19017 -5.62785 15
        commitDisps: 0 0 0
        Velocities   : 0 0 0
        commitAccels: 0 0 0
         unbalanced Load: 0 0 0
        ID : 0 1 2


ELEMENT DATA: NumEle: 200

Volume/Pressure Eight Node BbarBrick
Element Number: 0
Node 1 : 0
Node 2 : 1
Node 3 : 11
Node 4 : 10
Node 5 : 50
Node 6 : 51
Node 7 : 61
Node 8 : 60
Material Information :
 ElasticIsotropic3D
        tag: 0
        E: 0.3
        v: 4
        rho: 0
...
SP_Constraints: numConstraints: 150
SP_Constraint: 0         Node: 0 DOF: 0 value: 0
SP_Constraint: 1         Node: 0 DOF: 1 value: 0
SP_Constraint: 2         Node: 0 DOF: 2 value: 0
SP_Constraint: 3         Node: 1 DOF: 0 value: 0
SP_Constraint: 4         Node: 1 DOF: 1 value: 0
...
SP_Constraint: 149       Node: 49 DOF: 2 value: 0

MP_Constraints: numConstraints: 0

LOAD PATTERNS: num Patterns: 1

UniformExcitation - Need to Print the GroundMotion
Now I ask: what am I doing wrong? I mean, is the code for ground motion right? This is what I'm doing:

Code: Select all

theSeries = new TrigSeries (0, 200, 0.1, 0, 20);
GroundMotion *theMotion = new GroundMotion(0, 0, theSeries);
theLoadPattern = new UniformExcitation(*theMotion, 2, 1);
theDomain->addLoadPattern(theLoadPattern);
fmk
Site Admin
Posts: 5884
Joined: Fri Jun 11, 2004 2:33 pm
Location: UC Berkeley
Contact:

Post by fmk »

if you look at the element o/p you see that rho is still 0.0 .. IT SHOULD NOT BE .. now it will work once you get everything correct .. watch the period of sine versus period of structure .. here is the tcl script for a circular column with bbarBrick, ElasticIsotropic material, Trig series uniform excitation (which if i understand it is the model you are trying to do!) .. as boris said if you are just modelling something use the interpreter!

# Create ModelBuilder with 3 dimensions and 6 DOF/node
model basic -ndm 3 -ndf 3

# create the material
nDMaterial ElasticIsotropic 1 100 0.25 1.27

# Define geometry
# ---------------

# define some parameters
set eleArgs "1"

#set element stdBrick
set element bbarBrick

set nz 3
set nx 4
set ny 4

set nn [expr ($nz+1)*($nx+1)*($ny+1) ]

set root2 [expr sqrt(2.0)]
set -root2 [expr -$root2]

# mesh generation circular column root2 radius, height 10
block3D $nx $ny $nz 1 1 $element $eleArgs "
1 -1 -1 0
2 1 -1 0
3 1 1 0
4 -1 1 0
5 -1 -1 10
6 1 -1 10
7 1 1 10
8 -1 1 10
13 0 -$root2 0
14 $root2 0 0
15 0 $root2 0
16 -$root2 0 0
18 0 -$root2 10
19 $root2 0 10
20 0 $root2 10
21 -$root2 0 10
23 0 -$root2 5
24 $root2 0 5
25 0 $root2 5
26 -$root2 0 5
"


# set uniform load pattern
set lambda1 [eigen 1]
set period [expr 2.0 * 3.14159 / sqrt($lambda1)]
set tFinal [expr 10.0 * $period]

set accelSeries "Trig 0 $tFinal $period 0.0 -factor 0.2"
pattern UniformExcitation 1 1 -accel $accelSeries

# boundary conditions
fixZ 0.0 1 1 1

# -------------------------------------------------------------
# Create Some Recorders (graphics works linux, not sure windows)
# --------------------------------------------------------------

recorder Node -file Node.out -time -node $nn -dof 2 disp
recorder plot Node.out CenterNodeDisp 625 10 625 450 -columns 1 2
recorder display ShakingColumn 0 0 300 300 -wipe
prp -100 100 120.5
vup 0 1 0
display 1 0 1

# ---------------------------------------
# Create and Perform the dynamic analysis
# ---------------------------------------

# Create the transient analysis
test EnergyIncr 1.0e-10 20 0
algorithm Newton
numberer RCM
constraints Plain
integrator Newmark 0.5 0.25
analysis Transient

# Perform the transient analysis (20 sec)
# numSteps dt
set dt [expr $period/100.0]
set numIter [expr round($tFinal/$dt)]
analyze $numIter $dt

print ele 1
print node $nn
senseiwa
Posts: 8
Joined: Fri Nov 19, 2004 3:41 am
Location: Univ. Roma Tre

Post by senseiwa »

I want to develop in C++, as we want to integrate our software with the opensees FEM capabilities. The TCL interpreter is useless integrating C/C++ softwares. Now, I have some questions for translating your code into C++.

What does this mean?

Code: Select all

set lambda1 [eigen 1]
Is this true?

Code: Select all

pattern UniformExcitation 1 1 -accel $accelSeries 
    EQUALS TO
theSeries = new TrigSeries (0, 200, 2, 0, 0.2);
TrapezoidalTimeSeriesIntegrator *theTimeIntegrator = new TrapezoidalTimeSeriesIntegrator();
GroundMotion *theMotion = new GroundMotion(0, 0, accelSeries);
theMotion->setIntegrator(theTimeIntegrator);
theLoadPattern = new UniformExcitation(*theMotion, 1, 1);
theDomain->addLoadPattern(theLoadPattern);
Is this equivalent to use the C++ code I show for all the nodes whose Z coordinate is zero?

Code: Select all

fixZ 0.0 1 1 1 
    EQUALS TO \forall node \in Nodes
theConstraint = new SP_Constraint(constraint_id, node, 0, 0);
theConstraint = new SP_Constraint(constraint_id, node, 1, 0);
theConstraint = new SP_Constraint(constraint_id, node, 2, 0);
What about the time integrator? What does your ``model basic'' do?

I really need C++! :)
Post Reply