Skip to content

Commit

Permalink
Merge pull request #2360 from KratosMultiphysics/feature-poro-analysis
Browse files Browse the repository at this point in the history
Updating Poromechanics analysis and solver to new standard
  • Loading branch information
ipouplana authored Jun 23, 2018
2 parents 36c9fa2 + 4420597 commit c4c4fc3
Show file tree
Hide file tree
Showing 15 changed files with 1,114 additions and 627 deletions.
16 changes: 8 additions & 8 deletions applications/PoromechanicsApplication/README.md
Original file line number Diff line number Diff line change
Expand Up @@ -4,17 +4,17 @@ The Poromechanics Application contains developments in coupled solid-pore fluid

### Features:

- UPw small displacement element for saturated porous media (with
- UPw small displacement element for saturated porous media (with
equal order interpolation, unstable under incompressible-undrained
conditions)

- Stable UPw small displacement element for saturated porous media
- Stable UPw small displacement element for saturated porous media
(with higher order interpolation for displacements)

- FIC-Stabilized UPw small displacement element for saturated porous media
- FIC-Stabilized UPw small displacement element for saturated porous media
(with equal order interpolation for displacements)

- UPw Quasi-zero-thickness interface elements for defining cracks and
- UPw Quasi-zero-thickness interface elements for defining cracks and
joints

- Local linear elastic damage model (Simo-Ju and modified Von Mises)
Expand All @@ -24,8 +24,8 @@ Mises)

- Bilinear cohesive fracture model (for quasi-zero-thickness interface elements)

- Fracture propagation utility based on the combination of the
damage model with the insertion of interface elements after remeshing
- Fracture propagation utility based on the combination of the
damage model with the insertion of interface elements after remeshing
with GiD


Expand Down Expand Up @@ -61,8 +61,8 @@ Make sure that the following lines are properly set in the configure.sh file:
>
> -DUSE_PORO_MPI=ON \\
Uncomment the following line in
[poromechanics_main.py](https://github.com/KratosMultiphysics/Kratos/blob/master/applications/PoromechanicsApplication/custom_problemtype/Poromechanics_Application.gid/poromechanics_main.py) and [poromechanics_fracture_main.py](https://github.com/KratosMultiphysics/Kratos/blob/master/applications/PoromechanicsApplication/custom_problemtype/Poromechanics_Application.gid/poromechanics_fracture_main.py).
Uncomment the following line in
[KratosPoromechanics.py](https://github.com/KratosMultiphysics/Kratos/blob/master/applications/PoromechanicsApplication/custom_problemtype/Poromechanics_Application.gid/KratosPoromechanics.py).

> import KratosMultiphysics.TrilinosApplication as TrilinosApplication
Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,19 @@
from __future__ import print_function, absolute_import, division #makes KratosMultiphysics backward compatible with python 2.6 and 2.7

import KratosMultiphysics
import KratosMultiphysics.ExternalSolversApplication
# import KratosMultiphysics.TrilinosApplication as TrilinosApplication
import KratosMultiphysics.FluidDynamicsApplication
import KratosMultiphysics.SolidMechanicsApplication
import KratosMultiphysics.PoromechanicsApplication

from poromechanics_analysis import PoromechanicsAnalysis

if __name__ == "__main__":

with open("ProjectParameters.json",'r') as parameter_file:
parameters = KratosMultiphysics.Parameters(parameter_file.read())

model = KratosMultiphysics.Model()
simulation = PoromechanicsAnalysis(model,parameters)
simulation.Run()
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@ VALUE: 0.05
HELP: Length defining the domain of interaction between non-local damaged points.
QUESTION: Fracture_Propagation#CB#(true,false)
VALUE: false
HELP: Only for pre-defined crack-tips defined with interface elements.
HELP: Only for pre-defined crack tips defined with interface elements.
DEPENDENCIES: (true,RESTORE,Propagation_Length,#CURRENT#,RESTORE,Propagation_Damage,#CURRENT#,RESTORE,Propagation_Frequency,#CURRENT#,RESTORE,Propagation_Width,#CURRENT#,RESTORE,Propagation_Height,#CURRENT#,RESTORE,Correction_Tolerance,#CURRENT#)(false,HIDE,Propagation_Length,#CURRENT#,HIDE,Propagation_Damage,#CURRENT#,HIDE,Propagation_Frequency,#CURRENT#,HIDE,Propagation_Width,#CURRENT#,HIDE,Propagation_Height,#CURRENT#,HIDE,Correction_Tolerance,#CURRENT#)
QUESTION: Propagation_Damage
VALUE: 0.5
Expand All @@ -57,8 +57,8 @@ HELP: Number of steps after which the fracture propagation will be checked. Set


TITLE:Solver_Settings
QUESTION: Solution_Type#CB#(Quasi-Static,Dynamic)
VALUE: Quasi-Static
QUESTION: Solution_Type#CB#(quasi_static,Dynamic)
VALUE: quasi_static
QUESTION: Scheme_Type#CB#(Newmark)
VALUE: Newmark
QUESTION: Newmark_Beta
Expand All @@ -72,9 +72,9 @@ QUESTION: Rayleigh_Mass
VALUE: 0.0
QUESTION: Rayleigh_Stiffness
VALUE: 0.0
QUESTION: Strategy_Type#CB#(Newton-Raphson,Arc-Length)
VALUE: Newton-Raphson
DEPENDENCIES: (Newton-Raphson,HIDE,Desired_Iterations,#CURRENT#,HIDE,Max_Radius_Factor,#CURRENT#,HIDE,Min_Radius_Factor,#CURRENT#)(Arc-Length,RESTORE,Desired_Iterations,#CURRENT#,RESTORE,Max_Radius_Factor,#CURRENT#,RESTORE,Min_Radius_Factor,#CURRENT#)
QUESTION: Strategy_Type#CB#(newton_raphson,arc_length)
VALUE: newton_raphson
DEPENDENCIES: (newton_raphson,HIDE,Desired_Iterations,#CURRENT#,HIDE,Max_Radius_Factor,#CURRENT#,HIDE,Min_Radius_Factor,#CURRENT#)(arc_length,RESTORE,Desired_Iterations,#CURRENT#,RESTORE,Max_Radius_Factor,#CURRENT#,RESTORE,Min_Radius_Factor,#CURRENT#)
QUESTION: Convergence_Criterion#CB#(Displacement_criterion,Residual_criterion,And_criterion,Or_criterion)
VALUE: Displacement_criterion
DEPENDENCIES: (Displacement_criterion,RESTORE,Displacement_Relative_Tolerance,#CURRENT#,RESTORE,Displacement_Absolute_Tolerance,#CURRENT#,HIDE,Residual_Relative_Tolerance,#CURRENT#,HIDE,Residual_Absolute_Tolerance,#CURRENT#)(Residual_criterion,HIDE,Displacement_Relative_Tolerance,#CURRENT#,HIDE,Displacement_Absolute_Tolerance,#CURRENT#,RESTORE,Residual_Relative_Tolerance,#CURRENT#,RESTORE,Residual_Absolute_Tolerance,#CURRENT#)(And_criterion,RESTORE,Displacement_Relative_Tolerance,#CURRENT#,RESTORE,Displacement_Absolute_Tolerance,#CURRENT#,RESTORE,Residual_Relative_Tolerance,#CURRENT#,RESTORE,Residual_Absolute_Tolerance,#CURRENT#)(Or_criterion,RESTORE,Displacement_Relative_Tolerance,#CURRENT#,RESTORE,Displacement_Absolute_Tolerance,#CURRENT#,RESTORE,Residual_Relative_Tolerance,#CURRENT#,RESTORE,Residual_Absolute_Tolerance,#CURRENT#)
Expand Down
Original file line number Diff line number Diff line change
@@ -1,30 +1,30 @@
## GiD events --------------------------------------------------------------------------------------------------------------------------------------------------

proc InitGIDProject { dir } {

# Initialize ProblemType Menu
if { [GidUtils::IsTkDisabled] eq 0} {
if { [GidUtils::IsTkDisabled] eq 0} {
GiDMenu::Create "Poromechanics Application" PRE
GiDMenu::InsertOption "Poromechanics Application" [list "Parts"] 0 PRE "GidOpenConditions \"Parts\"" "" ""
GiDMenu::InsertOption "Poromechanics Application" [list "Dirichlet Constraints"] 1 PRE "GidOpenConditions \"Dirichlet_Constraints\"" "" ""
GiDMenu::InsertOption "Poromechanics Application" [list "Loads"] 2 PRE "GidOpenConditions \"Loads\"" "" ""
GiDMenu::InsertOption "Poromechanics Application" [list "Project Parameters"] 3 PRE "GidOpenProblemData" "" ""
GiDMenu::UpdateMenus
}

# Save ProblemTypePath
set ::Poromechanics_Application::ProblemTypePath $dir
}

#-------------------------------------------------------------------------------

proc AfterReadGIDProject { filename } {

# Save ProblemPath
set projectpath $filename
append projectpath .gid
set ::Poromechanics_Application::ProblemPath $projectpath

# Save ProblemName
# if {$::tcl_platform(platform) eq "windows"} {}
if {[regexp -all {\\} $filename] > 0} {
Expand All @@ -41,10 +41,10 @@ proc AfterReadGIDProject { filename } {
#-------------------------------------------------------------------------------

proc BeforeRunCalculation { batfilename basename dir problemtypedir gidexe args } {

# Set Parallel Configuration
set paralleltype [GiD_AccessValue get gendata Parallel_Configuration]

# Write Initial fractures data
if {([GiD_AccessValue get gendata Fracture_Propagation] eq true) && ($paralleltype ne "MPI")} {
# Define GiDPath
Expand All @@ -53,14 +53,14 @@ proc BeforeRunCalculation { batfilename basename dir problemtypedir gidexe args
regsub -all {\\} $gidexe {/} gidexe
}
set gidexe [string trimright $gidexe gid.exe]

if {[GiD_AccessValue get gendata Domain_Size] eq 2} {

source [file join $problemtypedir FracturePropagation2D.tcl]
WriteInitialFracturesData $dir $problemtypedir $gidexe

} else {

source [file join $problemtypedir FracturePropagation3D.tcl]
WriteInitialFracturesData $dir $problemtypedir $gidexe
}
Expand All @@ -73,14 +73,14 @@ proc BeforeRunCalculation { batfilename basename dir problemtypedir gidexe args
# Write ProjectParameters
source [file join $problemtypedir ProjectParameters.tcl]
WriteProjectParameters $basename $dir $problemtypedir $TableDict

# Copy python script in the problemdir
if {[GiD_AccessValue get gendata Fracture_Propagation] eq true} {
file copy -force [file join $problemtypedir poromechanics_fracture_main.py] [file join $dir MainKratos.py]
} else {
file copy -force [file join $problemtypedir poromechanics_main.py] [file join $dir MainKratos.py]
file copy -force [file join $problemtypedir KratosPoromechanics.py] [file join $dir MainKratos.py]
}

# Run the problem
set run 1
catch {
Expand All @@ -93,7 +93,7 @@ proc BeforeRunCalculation { batfilename basename dir problemtypedir gidexe args
Input files have been written.\n\
Run the case with: mpirun -np \[npartitions\] python3 MainKratos.py" ]]
}

### Measure time
#set start_time_1 [clock clicks]
#set end_time_1 [expr { [clock clicks]-$start_time_1 }]
Expand Down Expand Up @@ -128,7 +128,7 @@ proc Poromechanics_Application::PropagateFractures2D { } {
# Write new ProjectParameters
source [file join $::Poromechanics_Application::ProblemTypePath ProjectParameters.tcl]
WriteProjectParameters $::Poromechanics_Application::ProblemName $::Poromechanics_Application::ProblemPath $::Poromechanics_Application::ProblemTypePath $TableDict

# Quit GiD
GiD_Process Mescape Files Save
GiD_Process escape escape escape escape escape Quit
Expand All @@ -137,7 +137,7 @@ proc Poromechanics_Application::PropagateFractures2D { } {
#-------------------------------------------------------------------------------

proc Poromechanics_Application::PropagateFractures3D { } {

# Source Propagation Data and file
set PropagationData [source [file join $::Poromechanics_Application::ProblemPath PropagationData.tcl]]
source [file join $::Poromechanics_Application::ProblemTypePath FracturePropagation3D.tcl]
Expand All @@ -152,7 +152,7 @@ proc Poromechanics_Application::PropagateFractures3D { } {
# Write new ProjectParameters
source [file join $::Poromechanics_Application::ProblemTypePath ProjectParameters.tcl]
WriteProjectParameters $::Poromechanics_Application::ProblemName $::Poromechanics_Application::ProblemPath $::Poromechanics_Application::ProblemTypePath $TableDict

# Quit GiD
GiD_Process Mescape Files Save
GiD_Process escape escape escape escape escape Quit
Expand All @@ -161,60 +161,60 @@ proc Poromechanics_Application::PropagateFractures3D { } {
#--------------------------------------------------------------------------------------------------------------------------------------------------------------

proc Poromechanics_Application::CreateContactEntity { } {

set LayerName "Layer0"

## Contact Surface:

# set Line1 [list 111 1]
# set Line2 [list 110 0]

# # Orientation of Line1:
# # 0: SAME1ST. The normal to the line points towards the oppsite direction of the contact surface.
# # 1: DIFF1ST. The normal to the line points towards the contact surface.
# # Orientation of Line2:
# # 0: SAME1ST. The normal to the line points towards the contact surface.
# # 1: DIFF1ST. The normal to the line points towards the oppsite direction of the contact surface.

# GiD_Geometry create surface append contactsurface \
# $LayerName 2 $Line1 $Line2


## Contact Volume:

# set Surf1 [list 5 1]
# set Surf2 [list 186 0]

# # Orientation of Surf1:
# # 0: SAME1ST. The normal to the surface points towards the contact volume.
# # 1: DIFF1ST. The normal to the surface points towards the oppsite direction of the contact volume.
# # Orientation of Surf2:
# # 0: DIFF1ST. The normal to the surface points towards the oppsite direction of the contact volume.
# # 1: SAME1ST. The normal to the surface points towards the contact volume.

# Transformation matrix for a contact volume between Surf1 and Surf2. R means rotation, and T translation.

# set TransformMatrix [list Rxx Rxy Rxz Tx \
# Ryx Ryy Ryz Ty \
# Rzx Rzy Rzz Tz \
# 0.0 0.0 0.0 1.0]

# Transformation matrix for a contact volume between two equal surfaces, one over the other (zero-thickness interface elements)

# set TransformMatrix [list 1.0 0.0 0.0 0.0 \
# 0.0 1.0 0.0 0.0 \
# 0.0 0.0 1.0 0.0 \
# 0.0 0.0 0.0 1.0]

# Transformation matrix for a contact volume generated by a translation from Surf1 to Surf2

# set TransformMatrix [list 1.0 0.0 0.0 Tx \
# 0.0 1.0 0.0 Ty \
# 0.0 0.0 1.0 Tz \
# 0.0 0.0 0.0 1.0]

# Transformation matrix for a contact volume generated by a normal translation from Surf1 to Surf2

# # Vertex of Surf1
# set Vertex [GiD_Geometry get point 90]
# # Point1 of Surf1
Expand All @@ -226,7 +226,7 @@ proc Poromechanics_Application::CreateContactEntity { } {
# set TransformMatrix [NormalTranslationMatrix $Vertex $Point1 $Point2 $Distance]

# Transformation matrix for a contact volume generated by a rotation around a given axis from Surf1 to Surf2:

# # Point at the initial position of the rotation
# set InitPoint [GiD_Geometry get point 17]
# # Point at the final position of the rotation
Expand All @@ -238,7 +238,7 @@ proc Poromechanics_Application::CreateContactEntity { } {
# # Final point of the rotation Axis
# set FinalAxis [GiD_Geometry get point 19]
# set TransformMatrix [RotationMatrix $InitPoint $FinalPoint $Vertex $InitAxis $FinalAxis]

# GiD_Geometry create volume append $LayerName 2 \
# $Surf1 $Surf2 contactvolume $TransformMatrix
}
Expand All @@ -255,20 +255,20 @@ proc NormalTranslationMatrix {Vertex Point1 Point2 Distance} {
set Vy(0) [expr {[lindex $Point2 1]-[lindex $Vertex 1]}]
set Vy(1) [expr {[lindex $Point2 2]-[lindex $Vertex 2]}]
set Vy(2) [expr {[lindex $Point2 3]-[lindex $Vertex 3]}]
# Vector in local z direction (Cross product between Vx and Vy)

# Vector in local z direction (Cross product between Vx and Vy)
set Vz(0) [expr {$Vx(1)*$Vy(2)-$Vx(2)*$Vy(1)}]
set Vz(1) [expr {$Vx(2)*$Vy(0)-$Vx(0)*$Vy(2)}]
set Vz(2) [expr {$Vx(0)*$Vy(1)-$Vx(1)*$Vy(0)}]
set InvNorm [expr {1.0/sqrt($Vz(0)*$Vz(0)+$Vz(1)*$Vz(1)+$Vz(2)*$Vz(2))}]
set Vz(0) [expr {$Vz(0)*$InvNorm}]
set Vz(1) [expr {$Vz(1)*$InvNorm}]
set Vz(2) [expr {$Vz(2)*$InvNorm}]

set Tx [expr {$Distance*Vz(0)}]
set Ty [expr {$Distance*Vz(1)}]
set Tz [expr {$Distance*Vz(2)}]

return [list 1.0 0.0 0.0 $Tx \
0.0 1.0 0.0 $Ty \
0.0 0.0 1.0 $Tz \
Expand Down Expand Up @@ -303,7 +303,7 @@ proc RotationMatrix {InitPoint FinalPoint Vertex InitAxis FinalAxis} {
set A(0) [expr {$A(0)*$InvNorm}]
set A(1) [expr {$A(1)*$InvNorm}]
set A(2) [expr {$A(2)*$InvNorm}]

# Cosine of the angle of rotation
set CosAngle [expr {$Ri(0)*$Rf(0)+$Ri(1)*$Rf(1)+$Ri(2)*$Rf(2)}]
# Cross product between vectors Ri and Rf
Expand All @@ -312,7 +312,7 @@ proc RotationMatrix {InitPoint FinalPoint Vertex InitAxis FinalAxis} {
set n(2) [expr {$Ri(0)*$Rf(1)-$Ri(1)*$Rf(0)}]
# Sine of the angle of rotation (positive angle between 0º and 90º)
set SinAngle [expr {sqrt($n(0)*$n(0)+$n(1)*$n(1)+$n(2)*$n(2))}]

# Transformation Matrix
set Rxx [expr {$CosAngle+$A(0)*$A(0)*(1.0-$CosAngle)}]
set Rxy [expr {$A(0)*$A(1)*(1.0-$CosAngle)-$A(2)*$SinAngle}]
Expand All @@ -326,7 +326,7 @@ proc RotationMatrix {InitPoint FinalPoint Vertex InitAxis FinalAxis} {
set Tx [expr {[lindex $InitAxis 1]-($Rxx*[lindex $InitAxis 1]+$Rxy*[lindex $InitAxis 2]+$Rxz*[lindex $InitAxis 3])}]
set Ty [expr {[lindex $InitAxis 2]-($Ryx*[lindex $InitAxis 1]+$Ryy*[lindex $InitAxis 2]+$Ryz*[lindex $InitAxis 3])}]
set Tz [expr {[lindex $InitAxis 3]-($Rzx*[lindex $InitAxis 1]+$Rzy*[lindex $InitAxis 2]+$Rzz*[lindex $InitAxis 3])}]

return [list $Rxx $Rxy $Rxz $Tx \
$Ryx $Ryy $Ryz $Ty \
$Rzx $Rzy $Rzz $Tz \
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -16,15 +16,15 @@ export PYTHONPATH="$HOME/Kratos:$PYTHONPATH"

## Linux
# Setting PATHs for Kratos libs
export LD_LIBRARY_PATH="$HOME/Kratos/libs:$HOME/boost_1_62_0/stage/lib:$LD_LIBRARY_PATH"
export LD_LIBRARY_PATH="$HOME/Kratos/libs:$LD_LIBRARY_PATH"
# Setting PATH for Python 3
export PATH="/usr/bin/python3:$PATH"

## OS X
#~ # Setting PATHs for Kratos
#~ export DYLD_LIBRARY_PATH="$HOME/Kratos/libs:$HOME/boost_1_62_0/stage/lib:$DYLD_LIBRARY_PATH"
#~ # Setting PATH for Python 3
#~ export PATH="/Library/Frameworks/Python.framework/Versions/3.5/bin:${PATH}"
# # Setting PATHs for Kratos
# export DYLD_LIBRARY_PATH="$HOME/Kratos/libs:$DYLD_LIBRARY_PATH"
# # Setting PATH for Python 3
# export PATH="/Library/Frameworks/Python.framework/Versions/3.6/bin:${PATH}"

# Execute the program
"python3" "MainKratos.py" > "$1.info" 2> "$1.err"
Loading

0 comments on commit c4c4fc3

Please sign in to comment.