diff --git a/dG3D/README.txt b/dG3D/README.txt index 1bbb9c95a309f3863936245628987f5be3700883..89b4b36521a5927217786a965b6bfec0dd31e804 100644 --- a/dG3D/README.txt +++ b/dG3D/README.txt @@ -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 - diff --git a/dG3D/benchmarks/CMakeLists.txt b/dG3D/benchmarks/CMakeLists.txt index 893a03a833d710c51dc754d836a0a3e3bf4cb41f..07f97d664a40d4d5fdaff2e06e249a17c0be753c 100644 --- a/dG3D/benchmarks/CMakeLists.txt +++ b/dG3D/benchmarks/CMakeLists.txt @@ -264,3 +264,4 @@ add_subdirectory(multiScaleBodyForce_AnaPrev) add_subdirectory(honeycomb_compression) add_subdirectory(nonLinearTVE_uniaxial) add_subdirectory(nonLinearTVP_uniaxial) +add_subdirectory(nonLinearTVP_cube) diff --git a/dG3D/benchmarks/nonLinearTVP_cube/CMakeLists.txt b/dG3D/benchmarks/nonLinearTVP_cube/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..a0177858b739eed0c6a4bfbece73db3029b79186 --- /dev/null +++ b/dG3D/benchmarks/nonLinearTVP_cube/CMakeLists.txt @@ -0,0 +1,14 @@ +# test file + +set(PYFILE cubeTVP.py) + +set(FILES2DELETE + for*.csv + ene*.csv + Nod*.csv + disp* + stress* + IPVol* +) + +add_cm3python_test(${PYFILE} "${FILES2DELETE}") diff --git a/dG3D/benchmarks/nonLinearTVP_cube/cube.geo b/dG3D/benchmarks/nonLinearTVP_cube/cube.geo new file mode 100644 index 0000000000000000000000000000000000000000..e322959204440d14a93a5e085dce5e6ab743d194 --- /dev/null +++ b/dG3D/benchmarks/nonLinearTVP_cube/cube.geo @@ -0,0 +1,48 @@ +// 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) diff --git a/dG3D/benchmarks/nonLinearTVP_cube/cubeTVP.py b/dG3D/benchmarks/nonLinearTVP_cube/cubeTVP.py new file mode 100644 index 0000000000000000000000000000000000000000..d08bc850856e0655378442db67e8a6bd97189f23 --- /dev/null +++ b/dG3D/benchmarks/nonLinearTVP_cube/cubeTVP.py @@ -0,0 +1,200 @@ +#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) diff --git a/dG3D/benchmarks/nonLinearTVP_cube/relaxationSpectrumBetter_N30.csv b/dG3D/benchmarks/nonLinearTVP_cube/relaxationSpectrumBetter_N30.csv new file mode 100644 index 0000000000000000000000000000000000000000..ff29f4d9ddcee78818b14183cf907555f94e1208 --- /dev/null +++ b/dG3D/benchmarks/nonLinearTVP_cube/relaxationSpectrumBetter_N30.csv @@ -0,0 +1,61 @@ +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