units of vtk results

[ebrahimraeyat]

Hello.
As I asked here about modeling foundation on soil, I wanted to verify the results of xc. The results that getting from node information like getDisp is True, but when I want to view the results, I can't understand how can I set the correct units.

my input file:

# -*- coding: utf-8 -*-
from __future__ import print_function
'''Verification test taken from example Example-06 of 
   the CSI SAFE 2016 verification manual.'''


__author__= "Luis C. Pérez Tato (LCPT) and Ana Ortega (AOO)"
__copyright__= "Copyright 2015, LCPT and AOO"
__license__= "GPL"
__version__= "3.0"
__email__= "l.pereztato@gmail.com"

NumDivI= 8
NumDivJ= 8
CooMaxX= 30 * 12 # inch
CooMaxY= 20 * 12 # inch
E= 3000000 # Elastic modulus en lb/in2
nu= 0.2 # Poisson's ratio
G= 1250000
thickness= 15 # Cross section depth expressed in inches.
ptLoad= 100000 # Punctual load in lb.
soil_moduluse = 1000 # lb / in3
roz = 1 * soil_moduluse
display_results = True

import xc_base
import geom
import xc
from solution import predefined_solutions
from model import predefined_spaces
from materials import typical_materials
from model.boundary_cond.spring_bound_cond import ElasticFoundation


feProblem= xc.FEProblem()
preprocessor=  feProblem.getPreprocessor
nodes= preprocessor.getNodeHandler

# Problem type
modelSpace= predefined_spaces.StructuralMechanics3D(nodes)

# Define materials
memb1= typical_materials.defElasticMembranePlateSection(preprocessor, "memb1",E,nu,0.0,thickness)

seedElemHandler= preprocessor.getElementHandler.seedElemHandler
seedElemHandler.defaultMaterial= memb1.name
seedElemHandler.defaultTag= 1
elem= seedElemHandler.newElement("ShellMITC4",xc.ID([0,0,0,0]))

points= preprocessor.getMultiBlockTopology.getPoints
pt1= points.newPoint(geom.Pos3d(0.0,0.0,0.0))
pt2= points.newPoint(geom.Pos3d(CooMaxX,0.0,0.0))
pt3= points.newPoint(geom.Pos3d(CooMaxX,CooMaxY,0.0))
pt4= points.newPoint(geom.Pos3d(0.0,CooMaxY,0.0))
surfaces= preprocessor.getMultiBlockTopology.getSurfaces
surfaces.defaultTag= 1
s= surfaces.newQuadSurfacePts(pt1.tag,pt2.tag,pt3.tag,pt4.tag)
s.nDivI= NumDivI
s.nDivJ= NumDivJ
s.genMesh(xc.meshDir.I)

# Constraints
sides= s.getSides
# Edge iterator
for l in sides:
   if l.getEdge.getIVector.y == 0:
      for i in l.getEdge.getNodeTags():
         modelSpace.fixNode('00F_0FF', i)
   else:
      for i in l.getEdge.getNodeTags():
         modelSpace.fixNode('00F_F0F', i)


# Loads definition
lp0= modelSpace.newLoadPattern(name= 'DL')

node1= s.getNodeIJK(1, 1, 1)
node2= s.getNodeIJK(1, NumDivI+1, 1)
node3= s.getNodeIJK(1, 1, NumDivJ+1)
node4= s.getNodeIJK(1, NumDivI+1, NumDivJ+1)
nodes = (node1, node2, node3, node4)
for node in nodes:
  lp0.newNodalLoad(node.tag,xc.Vector([0,0,-ptLoad,0,0,0])) # Concentrated load


# nElems= s.getNumElements
# We add the load case to domain.
modelSpace.addLoadCaseToDomain(lp0.name)

foun = ElasticFoundation(soil_moduluse, roz, True)
setTotal= preprocessor.getSets["total"]
foun.generateSprings(setTotal) # set s as set did not work for display soil pressure

# Solution procedure
analysis= predefined_solutions.simple_static_linear(feProblem)
analOk= analysis.analyze(1)

node= s.getNodeIJK(1, int(NumDivI/2+1), int(NumDivJ/2+1))
print("Central node displacements: ", node.getDisp[2])
UZ= node1.getDisp[2]
print(f'Uz at corners is {UZ}')

if display_results:
   foun.displayPressures('soil Pressure', 1, 'lb / in2')
   from postprocess import output_handler
   from postprocess.config import default_config

   oh = output_handler.OutputHandler(modelSpace)
   # oh.getCameraParameters().viewName = "ZPos"
   # oh.displayFEMesh()
   # oh.displayBlocks()
   # oh.displayLocalAxes()
   # oh.displayLoads()
   # oh.displayReactions()
   # oh.displayDispRot(itemToDisp='uX')
   # oh.displayDispRot(itemToDisp='uY')
   oh.displayDispRot(itemToDisp='uZ')

Uz at corners is -0.05211809226666421 Central node displacements: 0.0007687480116149855

CSI SAFE and theoretical results is:
Uz at corners is -0.05308 Central node displacements: 0.00067

Also I don't know why my displacement result is not symmetric.

Thanks.

[lcpt]

Hi @ebrahimraeyat.

I don't understand the problem, you're using a 8x8 mesh (quite coarse) and you're getting 0.00077 in at center, and the theoretical displacement of the same point is 0.00067 in according to the table, so you get an error of 14.67%. SAFE gets a displacement of 0.00050 in for the same point and the same mesh, so the error is 25.37% (higher than yours).

For the corner point, you get a displacement of -0.052 in and the theoretical value is -0.05308, so the error in your model is 2.03%. SAFE gets a displacement of -0.06011 in which corresponds to an error of 13.24% (higher than yours again).

So, if I'm not wrong, you're using a 8x8 mesh and getting better results than SAFE, where is the problem?

With respect to the symmetry, for such a coarse mesh you are getting a nice symmetry there. If you use a finer mesh, you will see your results improve, but you will never obtain an absolutely perfect symmetry using a computer.

The units used for the graphics are controlled by the OutputUnits class here. If you want your results using imperial units, you'll need to modify the values on this object accordingly. We usually do something like this:

'''G1=graphical_reports.LoadCaseDispParameters(loadCaseName='GselfWeight',loadCaseDescr='G1: Poids propre',loadCaseExpr='1.0*GselfWeight',setsToDispLoads=[overallSet],setsToDispDspRot=[foundDeck,walls],setsToDispIntForc=[foundDeck,walls])
G1.unitsScaleLoads=1e-3
G1.unitsScaleDispl=1e3
G1.unitsDispl='[mm]'
G1.unitsScaleForc=1e-3
G1.unitsScaleMom=1e-3'''

to modify the factors that multiply displacements (and forces, ...) displayed with VTK.

[lcpt]

By the way, I've just seen the copyright of the code. Is this test really ours? I can't find it anywhere.

[ebrahimraeyat]

Thanks @lcpt for your response. I apologize for the delay in my reply. I’m not sure why I didn’t receive any email notifications about xcfem. Other repositories are okay.

Yes, it seems that the xc results are closer to the theoretical results than SAFE. However, my question was about the units of output. For example, the output shows -52.1 mm (mm or inch, it doesn’t matter), but it should display -0.00521. I think you divide the results by 1000 in your code. I need to investigate your method of setting units.

Regarding the symmetry of results, in SAFE it is really symmetrical. For this reason, I asked about it here to ensure that I didn’t make any mistakes in my code.

And regarding copyright, no, I just copied the entire code from a test and didn’t erase that part of the code.

[lcpt]

Hi, @ebrahimraeyat

About the units: by default, the postprocessor assumes that all the units in the analysis are the SI fundamental units (kg, m, s) and derived ones (N, m/s, m/s2, …). Then, to display the displacement in millimeters, it multiplies the displacements by 1000. You can check that the displacements returned by the nodes (using getDisp) are not affected.

Well, I don't know about SAFE, but I've seen the results of ANSYS for a lot of problems and I don't remember seeing a perfect symmetry ever. I think it's not important anyway.

If that's OK with you, you can change the copyright, and we can add this example to the set of verification tests.