Skip to content
Snippets Groups Projects
Commit 4e2c4cad authored by Ludovic Noels's avatar Ludovic Noels
Browse files

Merge branch 'ujwal3' into 'master' 

Ujwal3

See merge request !381
parents b7dc754d 7f6bf85a
No related branches found
No related tags found
1 merge request!381Ujwal3
Hide whitespace changes
Inline Side-by-side
dG3D/README.txt
+ 15
− 16
View file @ 4e2c4cad
......@@ -4,16 +4,16 @@ Build code from the command line
----------------------------------------------
In clusters, try to load CMake, OpenMPI, GCC, Python, SWIG
----------------------------------------------
----------------------------------------------
* create install directory
export CM3_DIR=$HOME/workspace
mkdir $CM3_DIR
* get PETSc and compile (mandatory if PETSc is not available)
* get PETSc and compile (mandatory if PETSc is not available)
#cloning source petsc
cd $CM3_DIR
cd $CM3_DIR
git clone -b main https://gitlab.com/petsc/petsc.git petsc
#choose PETSc version
......@@ -32,19 +32,19 @@ export PETSC_ARCH=arch-linux2-c-opt
make PETSC_DIR=$PETSC_DIR PETSC_ARCH=$PETSC_ARCH all
* code installation
* code installation
#get source code and gmsh
cd $CM3_DIR
git clone https://gitlab.onelab.info/gmsh/gmsh.git gmsh
# clone protected repository with your username
git clone https://username@gitlab.onelab.info/cm3/cm3Libraries.git cm3Libraries
git clone https://gitlab.onelab.info/gmsh/gmsh.git gmsh
# clone protected repository with your username
git clone https://username@gitlab.onelab.info/cm3/cm3Libraries.git cm3Libraries
#make a build directory
#make a build directory
cd $CM3_DIR/cm3Libraries/dG3D && mkdir release && cd release
# configuring
cmake -DCMAKE_BUILD_TYPE=Release \
# configuring
cmake -DCMAKE_BUILD_TYPE=Release \ (ccmake ..)
-DENABLE_MPI=ON \
-DENABLE_FLTK=OFF \
-DENABLE_EIGEN=OFF \
......@@ -54,7 +54,7 @@ cmake -DCMAKE_BUILD_TYPE=Release \
-DENABLE_WRAP_PYTHON=ON \
-DPETSC_DIR=$PETSC_DIR \
-DPETSC_ARCH=$PETSC_ARCH ..
# if there are conflicts with python library, include directory
# python include dir and library must be found
export PYTHON_INCLUDE_DIR=path-to-python-include
......@@ -74,7 +74,7 @@ cmake -DCMAKE_BUILD_TYPE=Release \
-DPETSC_DIR=$PETSC_DIR \
-DPETSC_ARCH=$PETSC_ARCH ..
# building code
make -j4
......@@ -83,9 +83,8 @@ export CM3_DIR=$HOME/workspace
export PATH=$CM3_DIR/cm3Libraries/dG3D/release/NonLinearSolver/gmsh:$PATH
export PYTHONPATH=$CM3_DIR/cm3Libraries/dG3D/release:$CM3_DIR/cm3Libraries/dG3D/release/NonLinearSolver/gmsh/utils/wrappers:$PYTHONPATH
# To use PETSc's built-in direct LU solver,
# To use PETSc's built-in direct LU solver,
# create the empty file .petscrc in $HOME and add the following options
-ksp_type preonly
-pc_type lu
-pc_factor_mat_solver_type petsc
dG3D/benchmarks/CMakeLists.txt
+ 1
− 0
View file @ 4e2c4cad
......@@ -264,3 +264,4 @@ add_subdirectory(multiScaleBodyForce_AnaPrev)
add_subdirectory(honeycomb_compression)
add_subdirectory(nonLinearTVE_uniaxial)
add_subdirectory(nonLinearTVP_uniaxial)
add_subdirectory(nonLinearTVP_cube)
dG3D/benchmarks/nonLinearTVP_cube/CMakeLists.txt 0 → 100644
+ 14
− 0
View file @ 4e2c4cad
# test file
set(PYFILE cubeTVP.py)
set(FILES2DELETE
for*.csv
ene*.csv
Nod*.csv
disp*
stress*
IPVol*
)
add_cm3python_test(${PYFILE} "${FILES2DELETE}")
dG3D/benchmarks/nonLinearTVP_cube/cube.geo 0 → 100644
+ 48
− 0
View file @ 4e2c4cad
// Cube.geo
mm=1.0e-3; // Units
n=1; // Number of layers (in thickness / along z)
m = 1; // Number of element + 1 along y
p = 1; // Number of element + 1 along x
L=1.*mm; // Cube lenght
sl1=L/n;
// Face z = 0 (point in anti-clockwise order)
Point(1)={0,0,0,sl1};
Point(2)={L,0,0,sl1};
Point(3)={L,L,0,sl1};
Point(4)={0,L,0,sl1};
Line(1)={1,2};
Line(2)={2,3};
Line(3)={3,4};
Line(4)={4,1};
Line Loop(5) = {3, 4, 1, 2};
Plane Surface(6) = {5};
// Cube extrusion
Extrude {0, 0, L} {
Surface{6}; Layers{n}; Recombine;
}
// Physical entities
// Volume
Physical Volume(29) = {1};
// Surface
Physical Surface(30) = {19}; // Left face x = 0
Physical Surface(31) = {27}; // Right face x = L
Physical Surface(32) = {6}; // Back face z = 0
Physical Surface(33) = {28}; // Front face z = L
Physical Surface(34) = {15}; // Down face y = 0
Physical Surface(35) = {23}; // Up face y = L
// Mesh
Transfinite Line {2, 4} = m Using Progression 1;
Transfinite Line {1, 3} = p Using Progression 1;
Transfinite Surface {6};
Recombine Surface{6};
// To save cube elongation
Physical Point(41) = {1}; // Point on left face (0,0,0)
Physical Point(42) = {4}; // Point on left face (0,L,0)
Physical Point(43) = {2}; // Point on right face (L,0,0)
Physical Point(44) = {3}; // Point on right face (L,L,0)
dG3D/benchmarks/nonLinearTVP_cube/cubeTVP.py 0 → 100644
+ 200
− 0
View file @ 4e2c4cad
#coding-Utf-8-*-
from gmshpy import *
from dG3Dpy import*
from math import*
import csv
import numpy as np
import os
# Load the csv data for relaxation spectrum
with open('relaxationSpectrumBetter_N30.csv', newline='') as csvfile:
reader = csv.reader(csvfile, delimiter=';')
relSpec = np.zeros((1,))
i = 0
for row in reader:
if i == 0:
relSpec[i] = float(' '.join(row))
else:
relSpec = np.append(relSpec, float(' '.join(row)))
i += 1
# print(relSpec)
N = int(0.5*(i-1))
relSpec[0:N+1] *= 1e-6 # convert to MPa
#script to launch beam problem with a python script
# material law
lawnum = 11 # unique number of law
E = relSpec[0] #MPa
nu = 0.3 #0.33
K = E/3./(1.-2.*nu) # Bulk mudulus
mu = E/2./(1.+nu) # Shear mudulus
rho = 7850e-9 #905e-12 # Bulk mass 905 kg/m3 -> 905e-12 tonne/mm3
Cp = rho*1900e+6 # 1900 J/kg/K -> 1900e+6 Nmm/tonne/K
KThCon = 0.14 # 0.14 W/m/K -> 0.14 Nmm/s/mm/K
Alpha = 0.6e-6 # 1/K
Tinitial = 273.15 + 70 # K
C1 = 136.17201827
C2 = 256.58268053 # K
# Temp Settings
Tref = 273.15+20.
TemFuncOpt = 0 # 0-> Shift Factor based derivatives of Ki; 1-> Gibson/Huang based derivatives
# Default tempfuncs
E_G = 2.80717267e+09 * 1e-6
E_R = 5.52466562e+08 * 1e-6
Cp_G = rho*1700e+6
Cp_R = rho*2100e+6
Alpha_G = 0.2e-6
Alpha_R = 1.e-6
KThCon_G = 0.10
KThCon_R = 1.5*0.12
k_gibson = 3.36435338e-02
m_gibson = 0.0
Tg = Tref
E_tempfunc = GlassTransitionScalarFunction(1, E_R/E_G, Tg, k_gibson, m_gibson, 0) # can't be used
Alpha_tempfunc = GlassTransitionScalarFunction(1, Alpha_R/Alpha_G, Tg, k_gibson, m_gibson, 0)
Cp_tempfunc = GlassTransitionScalarFunction(1, Cp_R/Cp_G, Tg, k_gibson, m_gibson, 0)
KThCon_tempfunc = GlassTransitionScalarFunction(1, KThCon_R/KThCon_G, Tg, k_gibson, m_gibson, 0)
ShiftFactor_tempfunc = negExpWLFshiftFactor(C1,C2,Tref)
sy0c = 10. #MPa
hc = 300.
sy0t = sy0c*0.78
ht = 300.
hardenc = LinearExponentialJ2IsotropicHardening(1, sy0c, hc, 0., 10.)
hardent = LinearExponentialJ2IsotropicHardening(2, sy0t, ht, 0., 10.)
hardenk = PolynomialKinematicHardening(3,1)
hardenk.setCoefficients(1,370)
law1 = NonLinearTVPDG3DMaterialLaw(lawnum,rho,E,nu,1e-6,Tinitial,Alpha,KThCon,Cp) # ,True,1e-8)
law1.setCompressionHardening(hardenc)
law1.setTractionHardening(hardent)
law1.setKinematicHardening(hardenk)
law1.setStrainOrder(11)
law1.setYieldPowerFactor(3.5)
law1.setNonAssociatedFlow(True)
law1.nonAssociatedFlowRuleFactor(0.5)
law1.setViscoelasticMethod(0)
law1.setTestOption(3)
law1.setShiftFactorConstantsWLF(C1,C2)
law1.setReferenceTemperature(Tref)
law1.setTempFuncOption(0)
law1.setViscoElasticNumberOfElement(N)
if N > 0:
for i in range(1, N+1):
law1.setViscoElasticData(i, relSpec[i], relSpec[i+N])
law1.setTemperatureFunction_ThermalDilationCoefficient(Alpha_tempfunc)
law1.setTemperatureFunction_ThermalConductivity(KThCon_tempfunc)
law1.setTemperatureFunction_SpecificHeat(Cp_tempfunc)
eta = constantViscosityLaw(1,3e4) #(1,3e4)
law1.setViscosityEffect(eta,0.15) # 0.21)
law1.setIsotropicHardeningCoefficients(1.,1.,1.)
# solver
sol = 2 # Gmm=0 (default) Taucs=1 PETsc=2
soltype = 1 # StaticLinear=0 (default) StaticNonLinear=1
nstep = 100 # number of step (used only if soltype=1)
ftime = 1 # Final time (used only if soltype=1)
tol=1.e-5 # relative tolerance for NR scheme (used only if soltype=1)
nstepArch=1 # Number of step between 2 archiving (used only if soltype=1)
fullDg = 0 #O = CG, 1 = DG
space1 = 0 # function space (Lagrange=0)
beta1 = 100
# geometry
geofile = "cube.geo" # "line.geo"
meshfile = "cube.msh" # "line.msh"
nfield = 29
pertLaw1 = dG3DMaterialLawWithTangentByPerturbation(law1,1e-8*(Tinitial)) # use with _pert*T0 (b/c In tangent_by_perturbation_TVP u use _pert*T0)
ThermoMechanicsEqRatio = 1.e1 # setConstitutiveExtraDofDiffusionEqRatio(ThermoMechanicsEqRatio); WHAT is this parameter?
thermalSource = True
mecaSource = True
myfield1 = ThermoMechanicsDG3DDomain(1000,nfield,space1,lawnum,fullDg,ThermoMechanicsEqRatio) # Added
myfield1.setConstitutiveExtraDofDiffusionAccountSource(thermalSource,mecaSource)
myfield1.stabilityParameters(beta1)
myfield1.ThermoMechanicsStabilityParameters(beta1,bool(1)) # Added
# creation of Solver
mysolver = nonLinearMechSolver(1000)
mysolver.createModel(geofile,meshfile,3,1)
mysolver.loadModel(meshfile)
mysolver.addDomain(myfield1)
mysolver.addMaterialLaw(law1)
# mysolver.addMaterialLaw(pertLaw1)
mysolver.Scheme(soltype)
mysolver.Solver(sol)
mysolver.snlData(nstep,ftime,tol)
mysolver.stepBetweenArchiving(nstepArch)
# BCs
umax = 0.001
fu_L = LinearFunctionTime(0,0,ftime,umax); # Linear Displacement with time
fu_L_2 = LinearFunctionTime(0,0,ftime,2*umax);
fu_L_3 = LinearFunctionTime(0,0,ftime,0.6*umax);
mysolver.displacementBC("Face",30,0,0.) # face x = 0 - Left Face fixed
mysolver.displacementBC("Face",31,0,fu_L) # face x = L - Right Face displaced - compression direction
mysolver.displacementBC("Face",31,1,fu_L_2)
mysolver.displacementBC("Face",31,1,fu_L_3)
mysolver.displacementBC("Face",34,1,0.) # face y = 0
mysolver.displacementBC("Face",32,2,0.) # face z = 0
# thermal BC
mysolver.initialBC("Volume","Position",nfield,3,Tinitial)
qmax = 0 #5.19e5 # Heat Flux in kN(ms)⁻¹(mm)⁻¹ [5.19e-1 W/m²]
fq_L = LinearFunctionTime(0,0,ftime,qmax); # Linear Flux with time
# mysolver.scalarFluxBC("Face",31,3,fq_L)
fT = LinearFunctionTime(0,Tinitial,ftime,70+Tinitial); # Linear Displacement with time
# fT = LinearFunctionTime(0,Tinitial,ftime,Tinitial); # Linear Displacement with time
mysolver.displacementBC("Volume",nfield,3,fT)
# Field-Output
mysolver.internalPointBuildView("strain_zz", IPField.STRAIN_ZZ, 1, 1);
mysolver.internalPointBuildView("svm",IPField.SVM, 1, 1)
mysolver.internalPointBuildView("sig_xx",IPField.SIG_XX, 1, 1)
mysolver.internalPointBuildView("sig_yy",IPField.SIG_YY, 1, 1)
mysolver.internalPointBuildView("sig_zz",IPField.SIG_ZZ, 1, 1)
mysolver.internalPointBuildView("sig_xy",IPField.SIG_XY, 1, 1)
mysolver.internalPointBuildView("sig_yz",IPField.SIG_YZ, 1, 1)
mysolver.internalPointBuildView("sig_xz",IPField.SIG_XZ, 1, 1)
mysolver.internalPointBuildView("epl",IPField.PLASTICSTRAIN, 1, 1)
mysolver.internalPointBuildView("temperature",IPField.TEMPERATURE, 1, 1)
mysolver.internalPointBuildView("qx",IPField.THERMALFLUX_X, 1, 1)
mysolver.internalPointBuildView("qy",IPField.THERMALFLUX_Y, 1, 1)
mysolver.internalPointBuildView("qz",IPField.THERMALFLUX_Z, 1, 1)
# Print/Archive History Outputs
# mysolver.archivingForceOnPhysicalGroup("Face", 31, 0, nstepArch) # Force on the right force
mysolver.archivingForceOnPhysicalGroup("Face", 30, 0, nstepArch) # Force on the left face === Reaction Force
# mysolver.archivingForceOnPhysicalGroup("Node", 41, 0, nstepArch) # Force on a point on the left face === Reaction Force
mysolver.archivingNodeDisplacement(43, 0, nstepArch) # Displacement of a point on the right face
mysolver.archivingIPOnPhysicalGroup("Volume",nfield, IPField.TEMPERATURE,IPField.MEAN_VALUE);
# mysolver.archivingForceOnPhysicalGroup("Face", 31, 3, nstepArch)
# mysolver.archivingIPOnPhysicalGroup("Face",31, IPField.THERMALFLUX_X,IPField.MEAN_VALUE);
# Solve
mysolver.solve()
check = TestCheck()
check.equal(-1.062888e-04,mysolver.getArchivedForceOnPhysicalGroup("Face", 30, 0),1.e-4)
dG3D/benchmarks/nonLinearTVP_cube/relaxationSpectrumBetter_N30.csv 0 → 100644
+ 61
− 0
View file @ 4e2c4cad
4.611542021776787639e+08
2.238302906718514562e+08
6.933795631409890950e+07
5.612462715236927569e+07
1.135300364704844803e+08
8.198232567065705359e+07
8.810335379910270870e+07
1.004036562814406157e+08
7.433283445541293919e+07
9.264592353174345195e+07
5.464365636767672747e+07
3.155253068493874371e+07
7.192803585835915804e+07
6.030332288085080683e+07
6.405368888230948895e+07
5.431639856282526255e+07
4.148370455185782164e+07
4.886423797939539701e+07
3.592593423976574838e+07
6.220685617107152194e+07
5.109722805920466781e+07
1.927525339279997349e+07
6.875568394042491913e+07
6.375761663832499832e+07
5.653656921271673590e+07
5.447126032106360048e+07
3.741191818174190074e+07
6.998646951026493311e+07
1.580391066094151512e+07
7.113430347675971687e+07
3.727696331510867923e+07
2.214070727389127921e-12
3.787529186956463470e-11
5.705652310353853072e-11
4.823563685188596367e-10
2.868546847643579377e-09
1.190414141551846151e-08
9.698182287690671800e-08
5.084635471550801825e-07
2.608335587038697655e-06
2.125642550198157774e-05
2.712603521128600229e-05
1.695535259534005365e-04
1.120501111743022913e-03
7.704679206158343736e-03
6.780508131648937953e-02
3.249279128365250568e-01
1.575706002496040981e+00
5.640665343564683631e+00
3.250919857455923534e+01
1.678403534861442665e+02
3.637830877100003022e+02
1.844418920184439685e+03
1.326507094467443676e+04
1.176262576267081022e+05
3.862514253418788430e+05
2.833541535623365082e+06
9.998967646092619747e+06
6.356490149844758213e+07
2.536265878395002186e+08
6.714104114903782606e+08
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment