Skip to content

Axisymmetric extension #1

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
migmolper wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Axisymmetric extension #1

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
a8549fe
First commit
Nov 26, 2022
13931a2
Second commit
Nov 27, 2022
b3590fa
Merge branch 'master' into axylsym
Dec 18, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
24 changes: 17 additions & 7 deletions nl-partsol/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,11 @@ ELSE ()
ENDIF ()

# Plain strain conditions
option(USE_PLAINSTRAIN "Use Plain-Strain conditions" ON)
option(USE_PLAIN_STRAIN "Use Plain-Strain conditions" ON)

# Axial symmetry conditions
option(USE_AXIAL_SYMMETRY "Use Axial-Symmetry conditions" OFF)


# Add LAPACK
option(USE_LAPACK "Use LAPACK solver library" ON)
Expand Down Expand Up @@ -51,11 +55,6 @@ link_libraries(m)
# Add PkgConfig
find_package(PkgConfig REQUIRED)

# Plain strain
if(USE_PLAINSTRAIN)
add_definitions("-DUSE_PLAINSTRAIN")
endif()

# MPI
if(USE_MPI)
find_package(MPI)
Expand Down Expand Up @@ -214,8 +213,19 @@ file(GLOB_RECURSE NL_PartSol_src

file(GLOB_RECURSE DATA ../build/*)

# Create executables
# Create executable (plain strain)
add_definitions(-DUSE_PLAIN_STRAIN=1 -DUSE_AXIAL_SYMMETRY=0 -DUSE_3D=0)
add_executable(${PROJECT_NAME} driver-nl-partsol.c ${NL_PartSol_src})
target_include_directories(${PROJECT_NAME} PUBLIC src)


# Create executable (Axial symmetry)
#add_definitions(-DUSE_PLAIN_STRAIN=0 -DUSE_AXIAL_SYMMETRY=1 -DUSE_3D=0)
#add_executable(${PROJECT_NAME}-axil driver-nl-partsol.c ${NL_PartSol_src})
#target_include_directories(${PROJECT_NAME}-axil PUBLIC src)


# Create executable (3D conditions)
#add_definitions(-DUSE_PLAIN_STRAIN=0 -DUSE_AXIAL_SYMMETRY=0 -DUSE_3D=1)
#add_executable(${PROJECT_NAME} driver-nl-partsol.c ${NL_PartSol_src})
#target_include_directories(${PROJECT_NAME} PUBLIC src)
2 changes: 2 additions & 0 deletions nl-partsol/driver-nl-partsol.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -68,6 +68,7 @@ bool Flag_Print_Convergence;
Load gravity_field;
bool Driver_EigenErosion;
bool Driver_EigenSoftening;
bool Driver_Fbar;
bool Petsc_Direct_solver;
bool Petsc_Iterative_solver;

Expand Down Expand Up @@ -97,6 +98,7 @@ int main(int argc, char *argv[]) {

// Default value for optional modulus
Driver_EigenErosion = false;
Driver_Fbar = false;

// Default values for the flags
Flag_Print_Convergence = false;
Expand Down
14 changes: 7 additions & 7 deletions nl-partsol/src/Constitutive/Constitutive.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@ int Stress_integration__Constitutive__(int p, Particle MPM_Mesh,
Input_SP.J = MPM_Mesh.Phi.Jbar.nV[p];
} else {
Input_SP.D_phi_n1 = MPM_Mesh.Phi.F_n1.nM[p];
Input_SP.J = MPM_Mesh.Phi.J_n1.nV[p];
Input_SP.J = MPM_Mesh.Phi.J_n1[p];
}

Output_SP = compute_Kirchhoff_Stress_Saint_Venant__Constitutive__(
Expand All @@ -46,7 +46,7 @@ int Stress_integration__Constitutive__(int p, Particle MPM_Mesh,
IO_State.Stress = MPM_Mesh.Phi.Stress.nM[p];
IO_State.W = &(MPM_Mesh.Phi.W[p]);
IO_State.D_phi_n1 = MPM_Mesh.Phi.F_n1.nM[p];
IO_State.J = MPM_Mesh.Phi.J_n1.nV[p];
IO_State.J = MPM_Mesh.Phi.J_n1[p];

STATUS =
compute_Kirchhoff_Stress_Hencky__Constitutive__(IO_State, MatProp_p);
Expand All @@ -68,7 +68,7 @@ int Stress_integration__Constitutive__(int p, Particle MPM_Mesh,
IO_State.J = MPM_Mesh.Phi.Jbar.nV[p];
} else {
IO_State.D_phi_n1 = MPM_Mesh.Phi.F_n1.nM[p];
IO_State.J = MPM_Mesh.Phi.J_n1.nV[p];
IO_State.J = MPM_Mesh.Phi.J_n1[p];
}

STATUS = compute_Kirchhoff_Stress_Neo_Hookean__Constitutive__(IO_State,
Expand All @@ -91,7 +91,7 @@ int Stress_integration__Constitutive__(int p, Particle MPM_Mesh,
IO_State.J = MPM_Mesh.Phi.Jbar.nV[p];
} else {
IO_State.D_phi_n1 = MPM_Mesh.Phi.F_n1.nM[p];
IO_State.J = MPM_Mesh.Phi.J_n1.nV[p];
IO_State.J = MPM_Mesh.Phi.J_n1[p];
}

// IO_State.Pressure = MPM_Mesh.Phi.lambda_pressure_n1.nV[p];
Expand Down Expand Up @@ -272,7 +272,7 @@ int stiffness_density__Constitutive__(int p, double *Stiffness_density,
// Compute local stiffness density
if (strcmp(MatProp_p.Type, "Neo-Hookean-Wriggers") == 0) {
IO_State.D_phi_n = MPM_Mesh.Phi.F_n.nM[p];
IO_State.J = MPM_Mesh.Phi.J_n1.nV[p];
IO_State.J = MPM_Mesh.Phi.J_n1[p];
STATUS = compute_stiffness_density_Neo_Hookean(
Stiffness_density, dN_alpha_n1, dN_beta_n1, dN_alpha_n, dN_beta_n,
IO_State, MatProp_p);
Expand All @@ -299,7 +299,7 @@ int stiffness_density__Constitutive__(int p, double *Stiffness_density,
IO_State.D_phi_n1 = MPM_Mesh.Phi.F_n1.nM[p];
IO_State.D_phi_n = MPM_Mesh.Phi.F_n.nM[p];
IO_State.dFdt = MPM_Mesh.Phi.dt_F_n1.nM[p];
IO_State.J = MPM_Mesh.Phi.J_n1.nV[p];
IO_State.J = MPM_Mesh.Phi.J_n1[p];
IO_State.alpha_4 = alpha_4;
STATUS = compute_stiffness_density_Newtonian_Fluid__Constitutive__(
Stiffness_density, dN_alpha_n1, dN_beta_n1, dN_alpha_n, dN_beta_n,
Expand Down Expand Up @@ -395,7 +395,7 @@ int compute_damage__Constitutive__(unsigned p, Particle MPM_Mesh,
const ChainPtr Beps_p = MPM_Mesh.Beps[p];
const double *kirchhoff_p = MPM_Mesh.Phi.Stress.nM[p];
const double *Strain_Energy_field = MPM_Mesh.Phi.W;
const double *J_n1 = MPM_Mesh.Phi.J_n1.nV;
const double *J_n1 = MPM_Mesh.Phi.J_n1;
const double *Vol_0 = MPM_Mesh.Phi.Vol_0.nV;

STATUS = Eigenerosion__Constitutive__(
Expand Down
16 changes: 6 additions & 10 deletions nl-partsol/src/Constitutive/Fluid/Newtonian-Fluid.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,8 +39,8 @@ int compute_Kirchhoff_Stress_Newtonian_Fluid__Constitutive__(

// Define and compute auxiliar tensor
#if NumberDimensions == 2
double L[4];
double E[4];
double L[5];
double E[5];
double Identity[4] = {1.0,0.0,0.0,1.0};
#else
double L[9];
Expand All @@ -55,11 +55,7 @@ int compute_Kirchhoff_Stress_Newtonian_Fluid__Constitutive__(
}
symmetrise__TensorLib__(E, L);

#if NumberDimensions == 2
double trace_E = E[0] + E[2];
#else
double trace_E = E[0] + E[4] + E[8];
#endif
double trace_E = I1__TensorLib__(E);

// Compute stress
for (unsigned i = 0; i < Ndim; i++) {
Expand Down Expand Up @@ -117,9 +113,9 @@ int compute_stiffness_density_Newtonian_Fluid__Constitutive__(

// Define auxiliar variables
#if NumberDimensions == 2
double L[4];
double E[4];
double b_n[4];
double L[5];
double E[5];
double b_n[5];
double E_dN_alpha_n1[2] = {0.0,0.0};
double E_dN_beta_n1[2] = {0.0,0.0};
double Lt_dN_alpha_n1[2] = {0.0,0.0};
Expand Down
2 changes: 1 addition & 1 deletion nl-partsol/src/Constitutive/Fluid/Newtonian-Fluid.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@
#include "Macros.h"
#include "Types.h"
#include "Matlib.h"
#include "Particles.h"
#include "Particles/compute-Strains.h"

/*!

Expand Down
6 changes: 3 additions & 3 deletions nl-partsol/src/Constitutive/Hyperelastic/Hencky.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,9 +43,9 @@ int compute_Kirchhoff_Stress_Hencky__Constitutive__(State_Parameters IO_State,
unsigned Ndim = NumberDimensions;

#if NumberDimensions == 2
double b[4];
double b[5];
double eigvec_b[4] = {0.0, 0.0, 0.0, 0.0};
double eigval_b[3] = {0.0, 0.0, 1.0};
double eigval_b[3] = {0.0, 0.0, 0.0};
double T_aux[4] = {0.0, 0.0, 0.0, 0.0};
#else
double b[9];
Expand Down Expand Up @@ -110,7 +110,7 @@ int compute_stiffness_density_Hencky__Constitutive__(
double G = E / (2.0 * (1.0 + nu));

#if NumberDimensions == 2
double b[4];
double b[5];
double eigvec_b[4] = {0.0, 0.0, 0.0, 0.0};
double eigval_b[2] = {0.0, 0.0};
double eigvec_T[4] = {0.0, 0.0, 0.0, 0.0};
Expand Down
2 changes: 1 addition & 1 deletion nl-partsol/src/Constitutive/Hyperelastic/Hencky.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@
#include "Macros.h"
#include "Types.h"
#include "Matlib.h"
#include "Particles.h"
#include "Particles/compute-Strains.h"
// clang-format on

/**
Expand Down
1 change: 1 addition & 0 deletions nl-partsol/src/Constitutive/Hyperelastic/Mooney-Rivlin.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
#include "Types.h"
#include "Matlib.h"
#include "Particles.h"
#include "Particles/compute-Strains.h"

/*******************************************************/

Expand Down
2 changes: 1 addition & 1 deletion nl-partsol/src/Constitutive/Hyperelastic/Neo-Hookean.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -109,7 +109,7 @@ int compute_stiffness_density_Neo_Hookean(

// Define auxiliar variables
#if NumberDimensions == 2
double b_n[4];
double b_n[5];
#else
double b_n[9];
#endif
Expand Down
2 changes: 1 addition & 1 deletion nl-partsol/src/Constitutive/Hyperelastic/Neo-Hookean.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@
#include "Macros.h"
#include "Types.h"
#include "Matlib.h"
#include "Particles.h"
#include "Particles/compute-Strains.h"
// clang-format on

/**
Expand Down
1 change: 1 addition & 0 deletions nl-partsol/src/Constitutive/Hyperelastic/Saint-Venant-Kirchhoff.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
#include "Types.h"
#include "Matlib.h"
#include "Particles.h"
#include "Particles/compute-Strains.h"

/*!

Expand Down
12 changes: 6 additions & 6 deletions nl-partsol/src/Formulations/Displacements-Pressure/U-p-Analisys.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -72,13 +72,13 @@ Fields allocate_Up_vars__Fields__(int NumParticles) {
#endif
}

Phi.J_n = allocZ__MatrixLib__(NumParticles, 1);
Phi.J_n = (double *)calloc(NumParticles,sizeof(double));

Phi.J_n1 = allocZ__MatrixLib__(NumParticles, 1);
Phi.J_n1 = (double *)calloc(NumParticles,sizeof(double));

for (int p = 0; p < NumParticles; p++) {
Phi.J_n.nV[p] = 1.0;
Phi.J_n1.nV[p] = 1.0;
Phi.J_n[p] = 1.0;
Phi.J_n1[p] = 1.0;
}

#if NumberDimensions == 2
Expand Down Expand Up @@ -169,8 +169,8 @@ void free_Up_vars__Fields__(Fields Phi) {
free__MatrixLib__(Phi.b_e_n1);
free__MatrixLib__(Phi.F_n);
free__MatrixLib__(Phi.F_n1);
free__MatrixLib__(Phi.J_n);
free__MatrixLib__(Phi.J_n1);
free(Phi.J_n);
free(Phi.J_n1);
free__MatrixLib__(Phi.dt_F_n);
free__MatrixLib__(Phi.dt_F_n1);
free__MatrixLib__(Phi.dt_DF);
Expand Down
6 changes: 3 additions & 3 deletions nl-partsol/src/Formulations/Displacements-Pressure/U-p-Newmark-beta.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -706,7 +706,7 @@ Define tributary nodes of the particle
/*
Compute Jacobian of the deformation gradient
*/
// MPM_Mesh.Phi.J_n1.nV[p] = I3__TensorLib__(F_n1_p);
// MPM_Mesh.Phi.J_n1[p] = I3__TensorLib__(F_n1_p);

/*
Free memory
Expand Down Expand Up @@ -1011,7 +1011,7 @@ compute_Volumetric_Constrain_Forces(Matrix Residual, Mask ActiveNodes,
/*
Get the Jacobian
*/
J_n1_p = MPM_Mesh.Phi.J_n1.nV[p];
J_n1_p = MPM_Mesh.Phi.J_n1[p];

/*
Impose the volumetric constrain
Expand Down Expand Up @@ -1690,7 +1690,7 @@ static void update_Particles(Nodal_Field D_up, Particle MPM_Mesh, Mesh FEM_Mesh,
/*
Replace the determinant of the deformation gradient
*/
MPM_Mesh.Phi.J_n.nV[p] = MPM_Mesh.Phi.J_n1.nV[p];
MPM_Mesh.Phi.J_n[p] = MPM_Mesh.Phi.J_n1[p];

/*
Iterate over the nodes of the particle
Expand Down
16 changes: 8 additions & 8 deletions nl-partsol/src/Formulations/Displacements-WaterPressure/U-pw-Analisys.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,13 +61,13 @@ Fields allocate_upw_vars__Fields__(int NumParticles) {
Allocate and initialise the Jacobian of the deformation gradient
ant its rate
*/
Phi.J_n = allocZ__MatrixLib__(NumParticles, 1);
Phi.J_n1 = allocZ__MatrixLib__(NumParticles, 1);
Phi.dJ_dt = allocZ__MatrixLib__(NumParticles, 1);
Phi.J_n = (double *)calloc(NumParticles,sizeof(double));
Phi.J_n1 = (double *)calloc(NumParticles,sizeof(double));
Phi.dJ_dt = (double *)calloc(NumParticles,sizeof(double));

for (int p = 0; p < NumParticles; p++) {
Phi.J_n.nV[p] = 1.0;
Phi.J_n1.nV[p] = 1.0;
Phi.J_n[p] = 1.0;
Phi.J_n1[p] = 1.0;
}

/*!
Expand Down Expand Up @@ -219,9 +219,9 @@ void free_upw_vars__Fields__(Fields Phi) {
free__MatrixLib__(Phi.dt_F_n);
free__MatrixLib__(Phi.dt_F_n1);
free__MatrixLib__(Phi.dt_DF);
free__MatrixLib__(Phi.J_n);
free__MatrixLib__(Phi.J_n1);
free__MatrixLib__(Phi.dJ_dt);
free(Phi.J_n);
free(Phi.J_n1);
free(Phi.dJ_dt);
free__MatrixLib__(Phi.Fbar);
free__MatrixLib__(Phi.Jbar);
free(Phi.W);
Expand Down
12 changes: 6 additions & 6 deletions nl-partsol/src/Formulations/Displacements-WaterPressure/U-pw-Newmark-beta.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -795,8 +795,8 @@ static void update_Local_State(
MPM_Mesh.Phi.rho_f.nV[p] *
MPM_Mesh.Phi.phi_f.nV[p]; /* Update density of the mixture */

MPM_Mesh.Phi.J_n1.nV[p] = J_n1_p; /* Update soil skeleton jacobian */
MPM_Mesh.Phi.dJ_dt.nV[p] =
MPM_Mesh.Phi.J_n1[p] = J_n1_p; /* Update soil skeleton jacobian */
MPM_Mesh.Phi.dJ_dt[p] =
dJ_dt_n1_p; /* Update soil skeleton rate of jacobian */

/*
Expand Down Expand Up @@ -1088,7 +1088,7 @@ static void compute_Rate_Mass_Fluid(Matrix Residual, Mask ActiveNodes,
/*
Get the rate of the jacobian for each material point
*/
dJ_dt_n1_p = MPM_Mesh.Phi.dJ_dt.nV[p];
dJ_dt_n1_p = MPM_Mesh.Phi.dJ_dt[p];

/*
Load intrinsic properties for the fluid phase to
Expand Down Expand Up @@ -1250,7 +1250,7 @@ static void compute_Flow_contribution_Fluid(Nodal_Field upw_n,
FT_n_p = transpose__TensorLib__(F_n_p);
Fm1_n1_p = Inverse__TensorLib__(F_n1_p);
FmT_n1_p = transpose__TensorLib__(Fm1_n1_p);
J_n1_p = MPM_Mesh.Phi.J_n1.nV[p];
J_n1_p = MPM_Mesh.Phi.J_n1[p];

/*
Get the reference volume for each material point (mixture)
Expand Down Expand Up @@ -1819,7 +1819,7 @@ static Matrix assemble_Tangent_Stiffness(Nodal_Field upw_n, Nodal_Field D_upw,
/*
Get the jacobian of the deformation gradient
*/
J_p = MPM_Mesh.Phi.J_n1.nV[p];
J_p = MPM_Mesh.Phi.J_n1[p];

/*
Get some usefull intermediate results
Expand Down Expand Up @@ -2386,7 +2386,7 @@ static void update_Particles(Nodal_Field D_upw, Particle MPM_Mesh,
/*
Replace the determinant of the deformation gradient
*/
MPM_Mesh.Phi.J_n.nV[p] = MPM_Mesh.Phi.J_n1.nV[p];
MPM_Mesh.Phi.J_n[p] = MPM_Mesh.Phi.J_n1[p];

/*
Compute the deformation energy (reference volume + material properties
Expand Down
Loading