Page MenuHomec4science
Files F107578686

No OneTemporary
Actions

File Metadata

Created
Wed, Apr 9, 03:11

View Options

diff --git a/docs/examples/theory/Dynamics_diagonalMass_1d/CMakeLists.txt b/docs/examples/theory/Dynamics_diagonalMass_1d/CMakeLists.txt
new file mode 100644
index 0000000..f114a59
--- /dev/null
+++ b/docs/examples/theory/Dynamics_diagonalMass_1d/CMakeLists.txt
@@ -0,0 +1,25 @@
+cmake_minimum_required(VERSION 3.1)
+
+# define a project name
+project(example)
+
+# set optimization level
+set(CMAKE_BUILD_TYPE Release)
+add_definitions(-DNDEBUG)
+
+# set C++ standard
+set(CMAKE_CXX_STANDARD 14)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+set(CMAKE_CXX_EXTENSIONS OFF)
+
+# load pkg-config
+find_package(PkgConfig)
+# load header-only modules using pkg-config
+pkg_check_modules(EIGEN3 REQUIRED eigen3)
+include_directories(${EIGEN3_INCLUDE_DIRS})
+
+# find pybind11
+find_package(pybind11 REQUIRED)
+
+# create the module
+pybind11_add_module(example example.cpp)
diff --git a/docs/examples/theory/Dynamics_diagonalMass_1d/damping.py b/docs/examples/theory/Dynamics_diagonalMass_1d/damping.py
new file mode 100644
index 0000000..508178f
--- /dev/null
+++ b/docs/examples/theory/Dynamics_diagonalMass_1d/damping.py
@@ -0,0 +1,79 @@
+
+import example
+import numpy as np
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+
+plt.style.use(['goose','goose-latex'])
+
+fig,ax = plt.subplots()
+
+# set variation in "nu" (critical damping at, nu == 1); colorbar settings
+exp = np.linspace(-3,+3,7)
+bnd = np.linspace(-3,+4,8)-.5
+lab = ['1/%d'%2**(-i) for i in exp if i<0] + ['%d'%2**i for i in exp if i>=0]
+nus = np.array([2**i for i in exp])
+cmap = plt.get_cmap('viridis',len(nus))
+
+# loop over all values of "nu"
+for inu,nu in enumerate(nus):
+
+ # set all parameters
+ nelem = 100
+ c = np.sqrt(nu)
+ G = 1.
+ rho = G / c**2.
+ eta = nu * rho
+ L = np.pi
+ h = L / nelem
+ qh = 2. * np.pi / h
+ qL = 2. * np.pi / L
+
+ # set external force: point force on the middle node
+ imid = int(nelem/2)
+ Fext = np.zeros((nelem+1)); Fext[imid] = -0.001 * G * h
+
+ # get time-scale : different from underdamped and overdamped systems
+ if ( nu/c < 2./qL ):
+ dt = 1 / ( c * qh )
+ dt /= 10.
+ else:
+ omega1 = ((nu*qh**2.)/2.) * np.sqrt( 2. - (2.*c/(qh*nu))**2. + np.sqrt(1.-(2.*c/(qh*nu))**2.) )
+ omega2 = ((nu*qh**2.)/2.) * np.sqrt( 2. - (2.*c/(qh*nu))**2. - np.sqrt(1.-(2.*c/(qh*nu))**2.) )
+ dt = max( 1./omega1 , 1./omega2 )
+
+ # set number of increments (total time constant for all samples)
+ T = 50.
+ ninc = int(T/dt)
+
+ # run simulation
+ u = example.compute(
+ rho = rho,
+ G = G,
+ eta = eta,
+ h = h,
+ Fext = Fext,
+ dt = dt,
+ ninc = ninc,
+ save_every = 1,
+ save_nodes = [imid],
+ )
+
+ # plot response
+ ax.plot(np.linspace(0,ninc*dt,len(u)),u/L,color=cmap(inu))
+
+# set plot labels
+plt.xlabel(r'$t$')
+plt.ylabel(r'$u$')
+
+# add colormap
+norm = mpl.colors.Normalize(vmin=min(exp),vmax=max(exp))
+sm = plt.cm.ScalarMappable(cmap=cmap,norm=norm)
+sm.set_array([])
+cbar = plt.colorbar(sm,ticks=exp,boundaries=bnd)
+cbar.set_ticklabels(lab)
+
+# show/save
+plt.show()
+plt.savefig('damping.svg')
+
diff --git a/docs/examples/theory/Dynamics_diagonalMass_1d/damping.svg b/docs/examples/theory/Dynamics_diagonalMass_1d/damping.svg
new file mode 100644
index 0000000..cd838f9
--- /dev/null
+++ b/docs/examples/theory/Dynamics_diagonalMass_1d/damping.svg
@@ -0,0 +1,21 @@
+<?xml version="1.0" encoding="utf-8" standalone="no"?>
+<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN"
+ "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
+<!-- Created with matplotlib (http://matplotlib.org/) -->
+<svg height="432pt" version="1.1" viewBox="0 0 576 432" width="576pt" xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink">
+ <defs>
+ <style type="text/css">
+*{stroke-linecap:butt;stroke-linejoin:round;}
+ </style>
+ </defs>
+ <g id="figure_1">
+ <g id="patch_1">
+ <path d="M 0 432
+L 576 432
+L 576 0
+L 0 0
+z
+" style="fill:none;"/>
+ </g>
+ </g>
+</svg>
diff --git a/docs/examples/theory/Dynamics_diagonalMass_1d/example.cpp b/docs/examples/theory/Dynamics_diagonalMass_1d/example.cpp
new file mode 100644
index 0000000..1f7ce15
--- /dev/null
+++ b/docs/examples/theory/Dynamics_diagonalMass_1d/example.cpp
@@ -0,0 +1,276 @@
+
+#define _USE_MATH_DEFINES // to use "M_PI" from "math.h"
+
+#include <iostream>
+#include <math.h>
+#include <Eigen/Eigen>
+#include <pybind11/pybind11.h>
+#include <pybind11/eigen.h>
+
+typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> MatD;
+typedef Eigen::Matrix<double, Eigen::Dynamic, 1, Eigen::ColMajor> ColD;
+typedef Eigen::Matrix<size_t, Eigen::Dynamic, 1, Eigen::ColMajor> ColS;
+
+// -------------------------------------------------------------------------------------------------
+
+ColD computeM(const ColD &x, double rho)
+{
+ ColD dNdxi(2), dNdx(2), xe(2), M(x.size());
+
+ // integration variables
+ double w = 1; // weight (N.B. such that total volume matches)
+ double J, V;
+
+ // shape function gradients : local coordinates
+ dNdxi(0) = -0.5;
+ dNdxi(1) = +0.5;
+
+ // extract number of elements : number of nodes - 1
+ size_t nelem = x.size() - 1;
+
+ // zero initialize force
+ M.setZero();
+
+ // loop over elements
+ for ( size_t e = 0 ; e < nelem ; ++e )
+ {
+ // element coordinates
+ xe(0) = x(e );
+ xe(1) = x(e+1);
+
+ // Jacobian
+ J = dNdxi(0) * xe(0) + dNdxi(1) * xe(1);
+
+ // element size
+ V = J * w;
+
+ // add to force
+ M(e ) += rho * V;
+ M(e+1) += rho * V;
+ }
+
+ return M;
+}
+
+// -------------------------------------------------------------------------------------------------
+
+ColD computeFu(const ColD &x, const ColD &u, double G)
+{
+ ColD dNdxi(2), dNdx(2), xe(2), ue(2), F(x.size());
+
+ // integration variables
+ double w = 2; // weight (N.B. xi == 0)
+ double J, V;
+
+ // local constitutive response
+ double sig, eps;
+
+ // shape function gradients : local coordinates
+ dNdxi(0) = -0.5;
+ dNdxi(1) = +0.5;
+
+ // extract number of elements : number of nodes - 1
+ size_t nelem = x.size() - 1;
+
+ // zero initialize force
+ F.setZero();
+
+ // loop over elements
+ for ( size_t e = 0 ; e < nelem ; ++e )
+ {
+ // element coordinates
+ xe(0) = x(e );
+ xe(1) = x(e+1);
+
+ // element displacements
+ ue(0) = u(e );
+ ue(1) = u(e+1);
+
+ // Jacobian
+ J = dNdxi(0) * xe(0) + dNdxi(1) * xe(1);
+
+ // element size
+ V = J * w;
+
+ // shape function gradients : global coordinates (constant over the element)
+ dNdx = ( 1./J ) * dNdxi;
+
+ // strain (constant over the element)
+ eps = dNdx(0) * ue(0) + dNdx(1) * ue(1);
+
+ // stress (constant over the element)
+ sig = G * eps;
+
+ // add to force
+ F(e ) += dNdx(0) * sig * V;
+ F(e+1) += dNdx(1) * sig * V;
+ }
+
+ return F;
+}
+
+// -------------------------------------------------------------------------------------------------
+
+ColD computeFv(const ColD &x, const ColD &v, double eta)
+{
+ ColD dNdxi(2), dNdx(2), xe(2), ve(2), F(x.size());
+
+ // integration variables
+ double w = 2; // weight (N.B. xi == 0)
+ double J, V;
+
+ // local constitutive response
+ double sig, epsdot;
+
+ // shape function gradients : local coordinates
+ dNdxi(0) = -0.5;
+ dNdxi(1) = +0.5;
+
+ // extract number of elements : number of nodes - 1
+ size_t nelem = x.size() - 1;
+
+ // zero initialize force
+ F.setZero();
+
+ // loop over elements
+ for ( size_t e = 0 ; e < nelem ; ++e )
+ {
+ // element coordinates
+ xe(0) = x(e );
+ xe(1) = x(e+1);
+
+ // element displacements
+ ve(0) = v(e );
+ ve(1) = v(e+1);
+
+ // Jacobian
+ J = dNdxi(0) * xe(0) + dNdxi(1) * xe(1);
+
+ // element size
+ V = J * w;
+
+ // shape function gradients : global coordinates (constant over the element)
+ dNdx = ( 1./J ) * dNdxi;
+
+ // strain (constant over the element)
+ epsdot = dNdx(0) * ve(0) + dNdx(1) * ve(1);
+
+ // stress (constant over the element)
+ sig = eta * epsdot;
+
+ // add to force
+ F(e ) += dNdx(0) * sig * V;
+ F(e+1) += dNdx(1) * sig * V;
+ }
+
+ return F;
+}
+
+// -------------------------------------------------------------------------------------------------
+
+MatD compute(
+ double rho, double G, double eta, double h, const ColD &Fext,
+ double dt, size_t ninc, size_t save_every, const ColS &save_nodes
+)
+{
+ size_t nel = static_cast<size_t>(Fext.size())-1;
+ double L = h * static_cast<double>(nel);
+ double t = 0.;
+ size_t isave = 0;
+ size_t jsave = static_cast<size_t>(std::floor(static_cast<double>(ninc)/static_cast<double>(save_every)));
+ size_t nsave = static_cast<size_t>(save_nodes.size());
+
+ MatD out(jsave,nsave);
+
+ ColD x (nel+1);
+ ColD u (nel+1);
+ ColD v (nel+1);
+ ColD v_n (nel+1);
+ ColD a (nel+1);
+ ColD a_n (nel+1);
+ ColD M (nel+1);
+ ColD Minv(nel+1);
+ ColD Fu (nel+1);
+ ColD Fv (nel+1);
+
+ x = ColD::LinSpaced( nel+1, 0.0, L );
+
+ u .setZero();
+ v .setZero();
+ a .setZero();
+ a_n.setZero();
+ v_n.setZero();
+
+ M = computeM(x,rho);
+ Minv = M.cwiseInverse();
+
+ for ( size_t inc = 0 ; inc < ninc ; ++inc )
+ {
+ // (1) new nodal positions (displacements)
+ // - apply update (nodal) : x_{n+1} = x_n + dt * v_n + .5 * dt^2 * a_n"
+ u += dt * v + ( .5 * std::pow(dt,2.) ) * a;
+ // - compute forces that are dependent on the displacement, but not on the velocity
+ Fu = computeFu(x,u,G);
+
+ // (2a) estimate nodal velocities
+ // - update velocities (DOFs)
+ v.noalias() = v_n + dt * a;
+ // - compute forces that are dependent on the velocity
+ Fv = computeFv(x,v,eta);
+ // - solve for accelerations (DOFs)
+ a.noalias() = Minv.cwiseProduct( Fext - Fu - Fv );
+ // - apply boundary conditions
+ a(0 ) = 0.;
+ a(nel) = 0.;
+ // - update velocities (DOFs)
+ v.noalias() = v_n + ( .5 * dt ) * ( a_n + a );
+ // - compute forces that are dependent on the velocity
+ Fv = computeFv(x,v,eta);
+
+ // (2b) new nodal velocities
+ // - solve for accelerations (DOFs)
+ a.noalias() = Minv.cwiseProduct( Fext - Fu - Fv );
+ // - apply boundary conditions
+ a(0 ) = 0.;
+ a(nel) = 0.;
+ // - update velocities (DOFs)
+ v.noalias() = v_n + ( .5 * dt ) * ( a_n + a );
+ // - compute forces that are dependent on the velocity
+ Fv = computeFv(x,v,eta);
+
+ // (3) new nodal accelerations
+ // - solve for accelerations (DOFs)
+ a.noalias() = Minv.cwiseProduct( Fext - Fu - Fv );
+ // - apply boundary conditions
+ a(0 ) = 0.;
+ a(nel) = 0.;
+
+ // store history
+ a_n = a; // accelerations (DOFs)
+ v_n = v; // nodal velocities
+ t += dt; // current time instance
+
+ // store output
+ if ( inc % save_every == 0 )
+ {
+ for ( size_t i=0; i<nsave; ++i ) out(isave,i) = u(save_nodes(i));
+ ++isave;
+ }
+ }
+
+ return out;
+}
+
+// -------------------------------------------------------------------------------------------------
+
+namespace py = pybind11;
+
+PYBIND11_MODULE(example,m)
+{
+ m.def(
+ "compute", &compute,
+ py::arg("rho"), py::arg("G"), py::arg("eta"), py::arg("h"), py::arg("Fext"), py::arg("dt"),
+ py::arg("ninc"), py::arg("save_every"), py::arg("save_nodes")
+ );
+}
+
diff --git a/docs/examples/theory/Dynamics_diagonalMass_1d/geometry.py b/docs/examples/theory/Dynamics_diagonalMass_1d/geometry.py
new file mode 100644
index 0000000..e24ce56
--- /dev/null
+++ b/docs/examples/theory/Dynamics_diagonalMass_1d/geometry.py
@@ -0,0 +1,62 @@
+
+import example
+import numpy as np
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+
+# set all parameters
+nelem = 100
+nu = 1./8.
+c = np.sqrt(nu)
+G = 1.
+rho = G / c**2.
+eta = nu * rho
+L = np.pi
+h = L / nelem
+qh = 2. * np.pi / h
+qL = 2. * np.pi / L
+
+# set external force: point force on the middle node
+imid = int(nelem/2)
+Fext = np.zeros((nelem+1)); Fext[imid] = -0.001 * G * h
+
+# get time-scale : different from underdamped and overdamped systems
+if ( nu/c < 2./qL ):
+ dt = 1 / ( c * qh )
+ dt /= 10.
+else:
+ omega1 = ((nu*qh**2.)/2.) * np.sqrt( 2. - (2.*c/(qh*nu))**2. + np.sqrt(1.-(2.*c/(qh*nu))**2.) )
+ omega2 = ((nu*qh**2.)/2.) * np.sqrt( 2. - (2.*c/(qh*nu))**2. - np.sqrt(1.-(2.*c/(qh*nu))**2.) )
+ dt = max( 1./omega1 , 1./omega2 )
+
+# set number of increments (total time constant for all samples)
+T = 50.
+ninc = int(T/dt)
+
+# run simulation
+u = example.compute(
+ rho = rho,
+ G = G,
+ eta = eta,
+ h = h,
+ Fext = Fext,
+ dt = dt,
+ ninc = ninc,
+ save_every = 1,
+ save_nodes = np.arange(nelem+1),
+)
+
+# plot response
+for i in range(u.shape[0])[::10]:
+
+ fig,ax = plt.subplots()
+ ax.plot(np.linspace(0,L,nelem+1),u[i,:]/L,color='k')
+
+ ax.axes.get_xaxis().set_visible(False)
+ ax.axes.get_yaxis().set_visible(False)
+
+ plt.xlim([0,L])
+ plt.ylim([-0.000012,0.000002])
+
+ plt.savefig('geometry_%05d.svg'%i)
+ plt.close()
diff --git a/docs/examples/theory/Dynamics_diagonalMass_1d/geometry_nu=0.125.mp4 b/docs/examples/theory/Dynamics_diagonalMass_1d/geometry_nu=0.125.mp4
new file mode 100644
index 0000000..cc29e05
Binary files /dev/null and b/docs/examples/theory/Dynamics_diagonalMass_1d/geometry_nu=0.125.mp4 differ

Event Timeline

c4science · Help