diff --git a/README.md b/README.md index 53e3961..07e46db 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # MPM-Py -A Material Point Method implementation using Python: MPM-Py. +A Material Point Method implementation using Python: MPM-Py. This program uses objected-oriented programming paradigm to represent and modeling the elements and its interaction in the material point method context. ## Installation @@ -10,47 +10,62 @@ git clone https://github.com/fabricix/MPM-Py.git ## Documentation -To create the documentation, run: +To create the documentation, install pdoc: ```bash -pdoc --html modules +pip3 install pdoc3 ``` - -To read the documentation, open +And then run ```bash -/html/modules/index.html +pdoc --html -c latex_math=True modules/ ``` +To read the documentation, open the file `/html/modules/index.html` using a web browser. + ## Requirements * Python 3.7.4 or superior * Matplotlib 3.3.4 or superior - ## Running tests and examples +In the folders `verification_problems` and `tests` there are examples showing and testing the functionalities of the program. + ### Mesh test +The file `tests/mesh-test.py` tests the mesh generations module by plotting the mesh and showing the number of elements, nodes and material points. + +Run this example as: + ```bash python mesh-test.py ``` ### Interpolation function test +The file `tests/interpolation_functions_test.py` shows the interpolation functions and its derivates over an one 1D element. + +Run this example as: ```bash python interpolation-functions-test.py ``` -### Single mass bar vibration problem +### Single mass vibration problem + +In this verification problem a single mass vibration is analyzed numerically and then the numerical solution is compared with the analytical one. +Run this example as: ```bash python mpm-single-mass-bar-vibration.py ``` ### Continuum bar vibration problem +In this verification problem a continuum bar vibration is analyzed numerically and then the numerical solution is compared with the analytical one. + +Run this example as: ```bash python mpm-continuum-bar-vibration.py -``` +``` \ No newline at end of file diff --git a/analytical/continuum_bar_vibration.py b/analitical_solutions/analitical_solution_continuum_bar_vibration.py similarity index 100% rename from analytical/continuum_bar_vibration.py rename to analitical_solutions/analitical_solution_continuum_bar_vibration.py diff --git a/analytical/single_mass_point_vibration.py b/analitical_solutions/analitical_solution_single_mass_vibration.py similarity index 100% rename from analytical/single_mass_point_vibration.py rename to analitical_solutions/analitical_solution_single_mass_vibration.py diff --git a/modules/element.py b/modules/element.py index 6d57e77..d2e0450 100644 --- a/modules/element.py +++ b/modules/element.py @@ -1,39 +1,38 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- -""" -Purpose -------- -Defines classes for finite elements representation -Data ----- -Created on Tue Mar 16 08:29:59 2021 +""" -Author ------- -Fabricio Fernandez, +Represent a finite element """ class bar_1D: - """ - Class to create 1D bar element with 2 nodes - - Attributes - ---------- - eid : int - Element identification - n1 : node type - Node 1 (left) - n2 : node type - Node 2 (right) - L : float - Element length - particles : list - List of particles in element - """ - def __init__(self): - self.eid = 0 # element id - self.n1 = 0 # node 1 - self.n2 = 0 # node 2 - self.L = 0 # element length - self.particles = [] # particles in element \ No newline at end of file + + """ + Represent an 1D finite element with 2 nodes + + Attributes + ---------- + + eid : int + Element identification + + n1 : node + Node 1 (left) + + n2 : node + Node 2 (right) + + L : float + Element length + + particles : list + List of particles in the element + """ + + def __init__(self): + self.eid = 0 # element id + self.n1 = 0 # node 1 (left) + self.n2 = 0 # node 2 (right) + self.L = 0 # element length + self.particles = [] # list of particles in element \ No newline at end of file diff --git a/modules/integration.py b/modules/integration.py index ca86ae7..e25da34 100644 --- a/modules/integration.py +++ b/modules/integration.py @@ -1,22 +1,21 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- + """ -Purpose -------- -Defines nodal integration functions -Data ----- -Created on Tue Mar 16 16:17:03 2021 +Defines nodal integration functions. +This functions in general have the objective to advance in time, or adding quantities at the nodes. -Author ------- -Fabricio Fernandez, """ def total_force_in_nodes(msh): """ Calculate total forces in nodes + + Arguments + --------- + msh: mesh + a mesh object """ for node in msh.nodes: node.f_tot = node.f_int + node.f_ext @@ -24,6 +23,14 @@ def total_force_in_nodes(msh): def momentum_in_nodes(msh,dt): """ Calculate momentum in nodes + + Arguments + --------- + msh: mesh + a mesh object + + dt: float + time step """ for inode in msh.nodes: inode.momentum += inode.f_tot*dt \ No newline at end of file diff --git a/modules/interpolation.py b/modules/interpolation.py index 46320bf..439bd8e 100644 --- a/modules/interpolation.py +++ b/modules/interpolation.py @@ -1,137 +1,160 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- -""" -Purpose -------- -Defines interpolation functions -Data ----- -Created on Mon Mar 15 15:25:44 2021 +r""" + +Defines interpolation functions in order to transfer quantities from particles to nodes. + +This functions, have the general form: + +.. math:: + f_I = \sum_p N_I(x_p) f_p + +where, +$$f_I$$: is a nodal quantity, +$$f_p$$: is a particle quantity, and +$$N_I(x_p)$$: is the weight or interpolation function of the node *I* evaluated at particle position $$x_p$$ -Author ------- -Fabricio Fernandez () """ def Ni(x,xI,L): - """ - Calculate the values of the interpolation function - - Arguments - --------- - x: float - point to calculate the function - xI: float - nodal point position - L: float - grid cell spacing - """ - - if(abs(x-xI)>=L): - return 0 - - if(-L<(x-xI) and (x-xI)<=0): - return 1+(x-xI)/L - - if(0<(x-xI) and (x-xI)=L): + return 0 + + if(-L<(x-xI) and (x-xI)<=0): + return 1+(x-xI)/L + + if(0<(x-xI) and (x-xI)=L): - return 0 - - if(-L<(x-xI) and (x-xI)<=0): - return 1/L - - if(0<(x-xI) and (x-xI)=L): + return 0 + + if(-L<(x-xI) and (x-xI)<=0): + return 1/L + + if(0<(x-xI) and (x-xI)) """ class linear_elastic: """ - Represents a linear elastic material. + Represents a linear elastic material Arguments --------- @@ -28,4 +21,4 @@ class linear_elastic: def __init__(self,E,density): self.E=E # Young's modulus - self.density=density # # density \ No newline at end of file + self.density=density # density \ No newline at end of file diff --git a/modules/mesh.py b/modules/mesh.py index 922e8b3..537c7d4 100644 --- a/modules/mesh.py +++ b/modules/mesh.py @@ -1,17 +1,10 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- + """ -Purpose -------- -Defines classes for finite element mesh representation -Data ----- -Created on Mon Mar 15 15:25:13 2021 +Represents a finite element mesh -Author ------- -Fabricio Fernandez () """ from modules import element @@ -82,7 +75,16 @@ def __init__(self,L,nelem): def put_particles_in_mesh(self,ppelem,material): """ - Function to distribute particles in elements mesh + Distributes particles in elements mesh + + Arguments + --------- + ppelem: int + number of particles per element + + material: material + a material object + """ self.ppelem=ppelem @@ -121,7 +123,12 @@ def put_particles_in_mesh(self,ppelem,material): def print_mesh(self,print_labels=True): """ - Function to print mesh in a plot + Function for print the mesh in a plot + + Arguments + --------- + print_labels: bool + determines if the label of the mesh will be plotted """ import matplotlib.pyplot as plt plt.cla() @@ -149,7 +156,7 @@ def print_mesh(self,print_labels=True): def print_mesh_info(self): """ - Function to print mesh informations + Function for print mesh informations """ print(20*'--') diff --git a/modules/node.py b/modules/node.py index 8a9c92d..f66cfd2 100644 --- a/modules/node.py +++ b/modules/node.py @@ -1,38 +1,39 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- + """ -Purpose -------- -Defines classes for finite element node representation -Data ----- -Created on Tue Mar 16 09:08:17 2021 +Represents a node in a finite element mesh -Author ------- -Fabricio Fernandez () """ + class node_1D: """ - Class to represents a 1D finite element node. + Represent a 1D node Attributes --------- nid : int node identification + x : float position + velocity : float nodal velocity + mass : float nodal mass + momentum : float nodal momentum (mass*velocity) + f_int : float nodal internal force + f_ext : float nodal external force + f_tot : float total force """ diff --git a/modules/particle.py b/modules/particle.py index d883940..1d0bd0b 100644 --- a/modules/particle.py +++ b/modules/particle.py @@ -1,22 +1,15 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- + """ -Purpose -------- -Defines classes for material point representation -Data ----- -Created on Mon Mar 15 15:22:44 2021 +Represents a material point -Author ------- -Fabricio Fernandez () """ class material_point: """ - Class to represents a material point. + Represent a material point. Attributes ---------- diff --git a/modules/solver.py b/modules/solver.py new file mode 100644 index 0000000..60ee2d9 --- /dev/null +++ b/modules/solver.py @@ -0,0 +1,147 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +r""" + +Perform the explicit solution in time of the motion equation: + +.. math:: + \frac{\partial \sigma_{ij}}{\partial x_j} + \rho b_i = \rho \ddot{u}_i + +where, +$$\sigma_{ij}$$: is the Cauchy stress tensor, +$$\rho$$: is the material density, and +$$b_i$$: is the body force per unit mass acting on the continuum, and +$$\ddot{u}_i$$: is the acceleration + +""" + +# local modules +from modules import interpolation as interpola # for interpolation tasks +from modules import integration as integra # for integration tasks +from modules import update # for updating tasks + +def explicit_solution(it,dt,time,msh,mpm_scheme,x_plot,y_plot,particle_plot,field_plot): + """ + Calculate the explicit solution of the motion equation using the MPM + + Arguments + --------- + it: float + Initial time + + dt: float + Time step + + time: float + total simulation time + + msh: mesh + a mesh object + + mpm_scheme: string + Update stress scheme ("USL", "MUSL" or "USF") + + x_plot: list + List to add the x variable of a result plot + + y_plot: list + List to add the y variable of a result plot + + particle_plot: int + Index of the particle to be plot the results + + field_plot: string + Field to be plotted ("velocity" or "position") + """ + + # main simulation loop + while it<=time: + + # update interpolation functions values + update.interpolation_functions_values(msh) + + # particle mass to grid nodal mass + interpola.mass_to_nodes(msh) + + # particle momentum to grid nodal momentum + interpola.momentum_to_nodes(msh) + + # impose essential boundary conditions (in fixed nodes set mv=0) + msh.elements[0].n1.momentum=0 + + # Update Stress First Scheme + if mpm_scheme=='USF': + + # calculate the grid nodal velocity + update.nodal_velocity(msh) + + # calculate particle strain increment + update.particle_strain_increment(msh,dt) + + # update particle density + update.particle_density(msh,dt) + + # update particle stress + update.particle_stress(msh,dt) + + # particle internal force to nodes + interpola.internal_force_to_nodes(msh) + + # particle external forces to nodes + interpola.external_force_to_nodes(msh) + + # calculate total force in node + integra.total_force_in_nodes(msh) + + # impose essential boundary conditions (in fixed nodes set f=m*a=0) + msh.elements[0].n1.f_tot=0 + + # integrate the grid nodal momentum equation + integra.momentum_in_nodes(msh,dt) + + # update particle velocity + update.particle_velocity(msh,dt) + + # update particle position + update.particle_position(msh,dt) + + # Modified Update Stress Last Scheme + if(mpm_scheme=='MUSL'): + + # recalculate the grid nodal momentum + update.nodal_momentum(msh) + + # impose essential boundary conditions (in fixed nodes v=0) + msh.elements[0].n1.velocity=0 + msh.elements[0].n1.momentum=0 + + # Modified Update Stress Last or Update Stress Last Scheme + if(mpm_scheme=='MUSL' or mpm_scheme=='USL'): + + # calculate the grid nodal velocity + update.nodal_velocity(msh) + + # calculate particle strain increment + update.particle_strain_increment(msh,dt) + + # update particle density + update.particle_density(msh,dt) + + # update particle stress + update.particle_stress(msh,dt) + + # reset all nodal values + update.reset_nodal_vaues(msh) + + # store data for plot + x_plot.append(it) + + if field_plot=='velocity': + y_plot.append(msh.particles[particle_plot].velocity) + + elif field_plot=='position': + y_plot.append(msh.particles[particle_plot].position) + + # advance in time + it+=dt diff --git a/modules/update.py b/modules/update.py index 4ad46a6..ca83383 100644 --- a/modules/update.py +++ b/modules/update.py @@ -1,23 +1,23 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- + """ -Purpose -------- -This file defines functions for updating tasks -Data ----- -Created on Tue Mar 16 16:24:45 2021 +This file defines functions for updating tasks -Author ------- -Fabricio Fernandez () """ from modules import interpolation as interp def particle_velocity(msh,dt): """ - update particle velocity + Update particle velocity + + Arguments + --------- + msh: mesh + a mesh object + dt: float + time step """ for ip in msh.particles: @@ -34,7 +34,14 @@ def particle_velocity(msh,dt): def particle_position(msh,dt): """ - update particle position + Update particle position + + Arguments + --------- + msh: mesh + a mesh object + dt: float + time step """ for ip in msh.particles: @@ -51,14 +58,24 @@ def particle_position(msh,dt): def nodal_velocity(msh): """ - calculate nodal velocity + Calculate nodal velocity + + Arguments + --------- + msh: mesh + a mesh object """ for inode in msh.nodes: inode.velocity=inode.momentum/inode.mass def nodal_momentum(msh): """ - calculate nodal momentum + Calculate nodal momentum + + Arguments + --------- + msh: mesh + a mesh object """ for inode in msh.nodes: inode.momentum=0 @@ -71,7 +88,14 @@ def nodal_momentum(msh): def particle_strain_increment(msh,dt): """ - calculate particle strain increment + Calculate particle strain increment + + Arguments + --------- + msh: mesh + a mesh object + dt: float + time step """ for ip in msh.particles: @@ -87,7 +111,15 @@ def particle_strain_increment(msh,dt): def particle_density(msh,dt): """ - update particle density + Update particle density + + Arguments + --------- + msh: mesh + a mesh object + dt: float + time step + """ for ip in msh.particles: @@ -95,7 +127,14 @@ def particle_density(msh,dt): def particle_stress(msh,dt): """ - update particle stress + Update particle stress + + Arguments + --------- + msh: mesh + a mesh object + dt: float + time step """ for ip in msh.particles: @@ -103,7 +142,12 @@ def particle_stress(msh,dt): def interpolation_functions_values(msh): """ - update the values of the nodal interpolation functions and its gradients + Update the values of the nodal interpolation functions and its gradients + + Arguments + --------- + msh: mesh + a mesh object """ for ip in msh.particles: @@ -120,7 +164,12 @@ def interpolation_functions_values(msh): def reset_nodal_vaues(msh): """ - reset all nodal values for the next step calculation + Reset all nodal values for the next step calculation + + Arguments + --------- + msh: mesh + a mesh object """ for inode in msh.nodes: inode.velocity = 0 diff --git a/mpm-continuum-bar-vibration.py b/mpm-continuum-bar-vibration.py deleted file mode 100755 index f1570f0..0000000 --- a/mpm-continuum-bar-vibration.py +++ /dev/null @@ -1,163 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Purpose -------- -This example approximates the continuum 1D bar vibration problem using MPM. - -Data ----- -Created on Sat Mar 20 20:25:53 2021 - -Author ------- -Fabricio Fernandez () -""" - -# external modules -import matplotlib.pyplot as plt # for plot -import numpy as np # for sin - -# local modules -from modules import mesh # for mesh definition -from modules import material # for material definition -from modules import interpolation as interpola # for interpolation tasks -from modules import integration as integra # for integration tasks -from modules import update # for updating tasks - -# bar length -L=25 - -# number of elements -nelements=15 - -# create an 1D mesh -msh = mesh.mesh_1D(L=L,nelem=nelements) - -# define a linear material -elastic = material.linear_elastic(E=100,density=1) - -# put particles in mesh element and set the material -msh.put_particles_in_mesh(ppelem=2,material=elastic) - -# simulation time -time = 60 # total time -dt = 0.1 # time step -it = 0.0 # initial time - -# verify time step -dt_critical=msh.elements[0].L/(elastic.E/elastic.density)**0.5 -dt = dt if dt < dt_critical else dt_critical - -# impose initial condition in particles -vo=0.1 -b1=np.pi/2.0/L -for ip in msh.particles: - ip.velocity=vo*np.sin(b1*ip.position) - -# variables for plot -x_plot=[] -y_plot=[] - -# MPM scheme integration -mpm_scheme='MUSL' # USL or USF or MUSL - -# main simulation loop -while it<=time: - - # update interpolation functions values - update.interpolation_functions_values(msh) - - # particle mass to grid nodal mass - interpola.mass_to_nodes(msh) - - # particle momentum to grid nodal momentum - interpola.momentum_to_nodes(msh) - - # impose essential boundary conditions (in fixed nodes set mv=0) - msh.elements[0].n1.momentum=0 - - # Update Stress First Scheme - if mpm_scheme=='USF': - - # calculate the grid nodal velocity - update.nodal_velocity(msh) - - # calculate particle strain increment - update.particle_strain_increment(msh,dt) - - # update particle density - update.particle_density(msh,dt) - - # update particle stress - update.particle_stress(msh,dt) - - # particle internal force to nodes - interpola.internal_force_to_nodes(msh) - - # particle external forces to nodes - interpola.external_force_to_nodes(msh) - - # calculate total force in node - integra.total_force_in_nodes(msh) - - # impose essential boundary conditions (in fixed nodes set f=m*a=0) - msh.elements[0].n1.f_tot=0 - - # integrate the grid nodal momentum equation - integra.momentum_in_nodes(msh,dt) - - # update particle velocity - update.particle_velocity(msh,dt) - - # update particle position - update.particle_position(msh,dt) - - # Modified Update Stress Last Scheme - if(mpm_scheme=='MUSL'): - - # recalculate the grid nodal momentum - update.nodal_momentum(msh) - - # impose essential boundary conditions (in fixed nodes v=0) - msh.elements[0].n1.velocity=0 - msh.elements[0].n1.momentum=0 - - # Modified Update Stress Last or Update Stress Last Scheme - if(mpm_scheme=='MUSL' or mpm_scheme=='USL'): - - # calculate the grid nodal velocity - update.nodal_velocity(msh) - - # calculate particle strain increment - update.particle_strain_increment(msh,dt) - - # update particle density - update.particle_density(msh,dt) - - # update particle stress - update.particle_stress(msh,dt) - - # reset all nodal values - update.reset_nodal_vaues(msh) - - # store for plot - x_plot.append(it) - y_plot.append(msh.particles[-1].velocity) - - # advance in time - it+=dt - -# plot mpm solution -plt.plot(x_plot,y_plot,'ob',markersize=2,label='mpm') - -# plot the analytical solution -from analytical import continuum_bar_vibration as cbv -[anal_xt,anal_vt, anal_t] = cbv.continuum_bar_vibration_solution(L,elastic.E,elastic.density,time,dt,vo,msh.particles[-1].position) -plt.plot(anal_t,anal_vt,'r',linewidth=2,label='analytical') - -# configure axis, legends and show plot -plt.xlabel('time (s)') -plt.ylabel('velocity (m/s)') -plt.legend() -plt.show() \ No newline at end of file diff --git a/mpm-single-mass-bar-vibration.py b/mpm-single-mass-bar-vibration.py deleted file mode 100755 index d7fe176..0000000 --- a/mpm-single-mass-bar-vibration.py +++ /dev/null @@ -1,160 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Purpose -------- -This example approximates the 1D single mass bar vibration problem using MPM. - -Data ----- -Created on Mon Mar 15 13:53:24 2021 - -Author ------- -Fabricio Fernandez () -""" - -# external modules -import matplotlib.pyplot as plt # for plot - -# local modules -from modules import mesh # for mesh definition -from modules import material # for material definition -from modules import interpolation as interpola # for interpolation tasks -from modules import integration as integra # for integration tasks -from modules import update # for updating tasks - -# bar length -L=1 - -# number of elements -nelements=1 - -# create an 1D mesh -msh = mesh.mesh_1D(L,nelements) - -# define a linear material -elastic = material.linear_elastic(E=50,density=1) - -# put particles in mesh element and set the material -msh.put_particles_in_mesh(ppelem=1,material=elastic) - -# simulation time -time = 10 # total time -dt = 0.001 # time step -it = 0 # initial time - -# verify time step -dt_critical=msh.elements[0].L/(elastic.E/elastic.density)**0.5 -dt = dt if dt < dt_critical else dt_critical - -# impose initial condition in particle -vo = 0.1 -msh.particles[-1].velocity=vo - -# variables for plot -x_plot=[] -y_plot=[] - -# MPM scheme integration -mpm_scheme='MUSL' # USL or USF or MUSL - -# main simulation loop -while it<=time: - - # update interpolation functions values - update.interpolation_functions_values(msh) - - # particle mass to grid nodal mass - interpola.mass_to_nodes(msh) - - # particle momentum to grid nodal momentum - interpola.momentum_to_nodes(msh) - - # impose essential boundary conditions (in fixed nodes set mv=0) - msh.elements[0].n1.momentum=0 - - # Update Stress First Scheme - if mpm_scheme=='USF': - - # calculate the grid nodal velocity - update.nodal_velocity(msh) - - # calculate particle strain increment - update.particle_strain_increment(msh,dt) - - # update particle density - update.particle_density(msh,dt) - - # update particle stress - update.particle_stress(msh,dt) - - # particle internal force to nodes - interpola.internal_force_to_nodes(msh) - - # particle external forces to nodes - interpola.external_force_to_nodes(msh) - - # calculate total force in node - integra.total_force_in_nodes(msh) - - # impose essential boundary conditions (in fixed nodes set f=m*a=0) - msh.elements[0].n1.f_tot=0 - - # integrate the grid nodal momentum equation - integra.momentum_in_nodes(msh,dt) - - # update particle velocity - update.particle_velocity(msh,dt) - - # update particle position - update.particle_position(msh,dt) - - # Modified Update Stress Last Scheme - if(mpm_scheme=='MUSL'): - - # recalculate the grid nodal momentum - update.nodal_momentum(msh) - - # impose essential boundary conditions (in fixed nodes v=0) - msh.elements[0].n1.velocity=0 - msh.elements[0].n1.momentum=0 - - # Modified Update Stress Last or Update Stress Last Scheme - if(mpm_scheme=='MUSL' or mpm_scheme=='USL'): - - # calculate the grid nodal velocity - update.nodal_velocity(msh) - - # calculate particle strain increment - update.particle_strain_increment(msh,dt) - - # update particle density - update.particle_density(msh,dt) - - # update particle stress - update.particle_stress(msh,dt) - - # reset all nodal values - update.reset_nodal_vaues(msh) - - # store for plot - x_plot.append(it) - y_plot.append(msh.particles[-1].position) - - # advance in time - it+=dt - -# plot mpm solution -plt.plot(x_plot,y_plot,'ob',markersize=2,label='mpm') - -# plot the analytical solution -from analytical import single_mass_point_vibration as smpv -[anal_xt, anal_t] = smpv.single_mass_point_vibration_solution(L,elastic.E,elastic.density,time,dt,L/2,vo) -plt.plot(anal_t,anal_xt,'r',linewidth=2,label='analytical') - -# configure axis, legends and show plot -plt.xlabel('time (s)') -plt.ylabel('displacement (m)') -plt.legend() -plt.show() \ No newline at end of file diff --git a/interpolation-functions-test.py b/tests/interpolation_functions_test.py similarity index 87% rename from interpolation-functions-test.py rename to tests/interpolation_functions_test.py index c5e5429..c55542b 100644 --- a/interpolation-functions-test.py +++ b/tests/interpolation_functions_test.py @@ -14,6 +14,10 @@ """ +# include the modules' path to the current path +import sys +sys.path.append("..") + # local modules from modules import mesh # for mesh generation from modules import interpolation # for interpolation diff --git a/mesh-test.py b/tests/mesh_test.py similarity index 88% rename from mesh-test.py rename to tests/mesh_test.py index 28f13f6..4f8b76c 100644 --- a/mesh-test.py +++ b/tests/mesh_test.py @@ -14,6 +14,10 @@ """ +# include the modules' path to the current path +import sys +sys.path.append("..") + # local modules from modules import mesh # for mesh generation from modules import material # for material definition diff --git a/verification_problems/mpm_continuum_bar_vibration.py b/verification_problems/mpm_continuum_bar_vibration.py new file mode 100755 index 0000000..9162007 --- /dev/null +++ b/verification_problems/mpm_continuum_bar_vibration.py @@ -0,0 +1,82 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Purpose +------- +This example approximates the continuum 1D bar vibration problem using MPM. + +Data +---- +Created on Sat Mar 20 20:25:53 2021 + +Author +------ +Fabricio Fernandez () +""" + +# include the modules' path to the current path +import sys +sys.path.append("..") + +# external modules +import matplotlib.pyplot as plt # for plot +import numpy as np # for sin + +# local modules +from modules import mesh # for mesh definition +from modules import material # for material definition +from modules import solver # for solving the problem in time + +# bar length +L=25 + +# number of elements +nelements=15 + +# create an 1D mesh +msh = mesh.mesh_1D(L=L,nelem=nelements) + +# define a linear material +elastic = material.linear_elastic(E=100,density=1) + +# put particles in mesh element and set the material +msh.put_particles_in_mesh(ppelem=2,material=elastic) + +# simulation time +time = 60 # total time +dt = 0.1 # time step +it = 0.0 # initial time + +# verify time step +dt_critical=msh.elements[0].L/(elastic.E/elastic.density)**0.5 +dt = dt if dt < dt_critical else dt_critical + +# impose initial condition in particles +vo=0.1 +b1=np.pi/2.0/L +for ip in msh.particles: + ip.velocity=vo*np.sin(b1*ip.position) + +# variables for plot +x_plot=[] +y_plot=[] + +# MPM scheme integration +mpm_scheme='MUSL' # USL or USF or MUSL + +# solve the problem in time +solver.explicit_solution(it,dt,time,msh,mpm_scheme,x_plot,y_plot,particle_plot=-1,field_plot='velocity') + +# plot mpm solution +plt.plot(x_plot,y_plot,'ob',markersize=2,label='mpm') + +# plot the analytical solution +from analitical_solutions import analitical_solution_continuum_bar_vibration as cbv +[anal_xt,anal_vt, anal_t] = cbv.continuum_bar_vibration_solution(L,elastic.E,elastic.density,time,dt,vo,msh.particles[-1].position) +plt.plot(anal_t,anal_vt,'r',linewidth=2,label='analytical') + +# configure axis, legends and show plot +plt.xlabel('time (s)') +plt.ylabel('velocity (m/s)') +plt.legend() +plt.show() \ No newline at end of file diff --git a/verification_problems/mpm_single_mass_vibration.py b/verification_problems/mpm_single_mass_vibration.py new file mode 100755 index 0000000..d0e56b8 --- /dev/null +++ b/verification_problems/mpm_single_mass_vibration.py @@ -0,0 +1,79 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Purpose +------- +This example approximates the 1D single mass bar vibration problem using MPM. + +Data +---- +Created on Mon Mar 15 13:53:24 2021 + +Author +------ +Fabricio Fernandez () +""" + +# include the modules' path to the current path +import sys +sys.path.append("..") + +# external modules +import matplotlib.pyplot as plt # for plot + +# local modules +from modules import mesh # for mesh definition +from modules import material # for material definition +from modules import solver # for solving the problem in time + +# bar length +L=1 + +# number of elements +nelements=1 + +# create an 1D mesh +msh = mesh.mesh_1D(L,nelements) + +# define a linear material +elastic = material.linear_elastic(E=50,density=1) + +# put particles in mesh element and set the material +msh.put_particles_in_mesh(ppelem=1,material=elastic) + +# simulation time +time = 10 # total time +dt = 0.001 # time step +it = 0 # initial time + +# verify time step +dt_critical=msh.elements[0].L/(elastic.E/elastic.density)**0.5 +dt = dt if dt < dt_critical else dt_critical + +# impose initial condition in particle +vo = 0.1 +msh.particles[-1].velocity=vo + +# variables for plot +x_plot=[] +y_plot=[] + +# MPM scheme integration +mpm_scheme='MUSL' # USL or USF or MUSL + +# solve the problem in time +solver.explicit_solution(it,dt,time,msh,mpm_scheme,x_plot,y_plot,particle_plot=-1,field_plot='position') + +# plot mpm solution +plt.plot(x_plot,y_plot,'ob',markersize=2,label='mpm') + +# plot the analytical solution +from analitical_solutions import analitical_solution_single_mass_vibration as smpv +[anal_xt, anal_t] = smpv.single_mass_point_vibration_solution(L,elastic.E,elastic.density,time,dt,L/2,vo) +plt.plot(anal_t,anal_xt,'r',linewidth=2,label='analytical') + +# configure axis, legends and show plot +plt.xlabel('time (s)') +plt.ylabel('displacement (m)') +plt.legend() +plt.show() \ No newline at end of file