Truss

examples/ex12/ex12.cpp

@@ -47,7 +47,12 @@ int main(int argc, char **argv) {
   // double beta = 0;
   // double gamma = 0.5;
 
+
   ShapeFunctions();
+
+  //printf("execut until here\n");
+  //exit (EXIT_FAILURE);
+
   ReadMaterialProperties();
   /*  Step-1: Calculate the mass matrix similar to that of belytschko. */
   Assembly((char *)"mass"); // Add Direct-lumped as an option

femtech_config.sh

@@ -0,0 +1,24 @@
+#!/bin/bash
+#
+#UNAMEX="ubuntu"
+#cd /home/$UNAMEX
+#sudo apt-get update
+#sudo apt-get install -y build-essential
+#sudo apt-get install -y cmake-curses-gui
+#sudo apt-get install -y libblas-dev liblapack-dev
+## sudo apt-get install -y openmpi-bin openmpi-common openssh-client openssh-server libopenmpi1.10
+#sudo apt-get install -y openmpi-bin openmpi-common libopenmpi-dev
+#git clone https://github.com/PSUCompBio/FemTech
+#cd FemTech
+mkdir build
+cd build
+cmake .. -DEXAMPLES=ON -DEXAMPLE12=ON -DEXAMPLE11=ON
+make -j8
+
+#make -j8 ex12 In the examples/ex12/ there is also a makefile. use it
+
+#mpirun -n 1 ex11 1-elt-cube.k > debug_log.txt
+#mpirun -n 1 ex12 1-elt-truss.k > debug_log.txt
+#paraview 1-elt-cube.pvd  here the *.pvd file
+
+#fileellelkfalfhf

include/FemTech.h

@@ -32,22 +32,29 @@ int LineToArray(const bool IntOrFloat, const bool CheckLastVal, \
     const char *Delim = " \t", void **Array = NULL);
 bool ReadInputFile(const char *FileName);
 bool PartitionMesh();
+
 void GaussQuadrature3D(int element, int nGaussPoint, double *Chi,double *GaussWeights);
+
 void ShapeFunctions();
 void ShapeFunction_C3D8(int e, int gp, double *Chi, double *detJ);
 void ShapeFunction_C3D4(int e, int gp, double *Chi, double *detJ);
 void ShapeFunction_T3D2(int e, int gp, double *Chi, double *detJ);
+void get_cos(int e, double *cx, double *cy, double *cz);
+
 void ReadMaterialProperties();
 void ReadBoundaryCondition(void);
+
 void AllocateArrays();
 
 void Assembly(char *operation);
 void StiffnessElementMatrix(double* Ke, int e);
 void MassElementMatrix(double* Me, int e);
+
 void WriteVTU(const char* FileName, int step, double time);
 void WritePVD(const char* FileName, int step, double time);
+
 void FreeArrays();
-void ReadMaterialProperties();
+
 void ApplySteadyBoundaryConditions(void);
 void SolveSteadyImplicit(void);
 void SolveUnsteadyNewmarkImplicit(double beta, double gamma, double dt, \

python_code/Main.py

@@ -0,0 +1,225 @@
+# ****************************************************************************
+# Name           : Main.py
+# Author         : Andres Valdez
+# Version        : 1.0
+# Description    : 3D truss element, pull-back matrix operations
+#                  This set of functions receives Nodal displaments, and
+#                  provides the transformation matrix, from global to local
+# Data	     : 05-05-2019
+# ****************************************************************************
+
+from __future__ import unicode_literals
+import os, sys
+import numpy as np
+import scipy as scp
+from scipy import optimize
+from scipy.linalg import solve
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+import matplotlib.mlab as mlab
+from math import pi, sin, sqrt, pow
+
+def get_cos(x,y,z,l0):
+    '''
+    Receives the variations of lengths between two consecutive mechanical
+    states. And computes the direction towards the initial state.
+    '''
+    cx = (x[1] - x[0]) / l0
+    cy = (y[1] - y[0]) / l0
+    cz = (z[1] - z[0]) / l0
+    return cx, cy, cz
+
+def map_global_to_local(cx,cy,cz):
+    '''
+    Here we define the 3d to 2d mapping. (x,y,z)--->(xo,xf)
+    '''
+    T = np.zeros((2,6))
+    T[0,0], T[0,1], T[0,2] = cx, cy, cz
+    T[1,3], T[1,4], T[1,5] = cx, cy, cz
+
+    return T
+
+def vec_global_to_local(T,f):
+    '''
+    Receives any global vector (elemental) and the transformation matrix.
+    Returns the local vector for Isoparametric implementation
+    '''
+    return np.matmul(T,f)
+
+def vec_local_to_global(T,f):
+    '''
+    Receives the local vector and the global to local mapping, and returns
+    the global vector.
+    '''
+    return np.matmul(T.T,f)
+
+def gauss_point(npts):
+    '''
+    This function returns the localization of a gauss points and weights
+    consideering the ref domain (-1,1).
+    '''
+    pg, wg = np.zeros(npts) , np.zeros(npts)
+
+    if(npts==1):
+        pg[0] = 0.0
+        wg[0] = 2
+    if(npts==2):
+        pg[0] , pg[1] = 1/sqrt(3) , -1/sqrt(3)
+        wg[0] , wg[1] = 1 , 1
+    return pg, wg
+
+def shape_func_local(npts):
+    '''
+    Here we define the local shape functions. Isoparametric (-1,1).
+    Just for a 2-node line
+    shl(0,i,l) = x-derivative of shape function i
+    shl(1,i,l) =  shape function defined at node i
+    l: gauss point
+    '''
+    shl = np.zeros((2,2,npts))
+
+    pg , wg = gauss_point(npts)
+
+    for k in range(npts):
+        shl[0,0,k] = -0.5
+        shl[1,0,k] = -0.5 * (pg[k] - 1.0)
+
+        shl[0,1,k] = 0.5
+        shl[1,1,k] = -0.5 * (pg[k] + 1.0)
+
+    return shl
+
+def shape_func_global(npts,x):
+    '''
+    Here we define the global shape functions
+    shg(0,i,l) = x-derivative of shape function i
+    shg(1,i,l) =  shape function defined at node i
+    x[0],x[1] local-coordinates of line nodes
+
+    '''
+    det = np.zeros(npts)
+    shg = np.zeros((2,2,npts))
+
+    shl = shape_func_local(npts)
+
+    for k in range(npts):
+        det[k] = sqrt((x[1]-x[0])**2)
+
+        for i in range(2):
+            shg[0,i,k] = shl[0,i,k] / det[k]
+            shg[1,i,k] = shl[1,i,k]
+
+    return shg,det
+
+def pos_truss(shg,x,y,z,l):
+    '''
+    Here we define the x,y,z-variables for integration @ each gauss-point
+    '''
+    px = 0.0
+
+    for k in range(2):
+        px = px + shg[1,k,l] * x[k]
+    return px,py,pz
+
+def elemental_force(force,x,npts):
+    '''
+    Here we define the elemental force vector.
+    Receives local forces, local coordinates
+    '''
+    fe = np.zeros(2)
+    pg , wg = gauss_point(npts)
+    shg , det = shape_func_global(npts,x)
+    for l in range(npts):
+        for k in range(2):
+            fe[k] = fe[k] + force[k] * wg[l] * shg[1,k,l] * det[l]
+
+    return fe
+
+def Green_Lagrange(npts,x,u):
+    '''
+    Here we define the elemental Green-Lagrange tensor.
+    Receives number of gauss points and local coordinates and a local
+    displacement.
+    '''
+    du = np.zeros(2)
+    Eu = np.zeros((2,2))
+    pg , wg = gauss_point(npts)
+    shg , det = shape_func_global(npts,x)
+    for l in range(npts):
+        for k in range(2):
+            du[k] = du[k] + shg[0,k,l] * wg[l] * det[l] * u[k]
+
+    Eu[0,0] = du[0] + 0.5 * du[0]**2
+    Eu[1,1] = du[1] + 0.5 * du[1]**2
+    return Eu
+
+def elemental_int_force(Eu,E,A):
+    '''
+    Receives the Green-Lagrange strain tensor, the Young modulus, and the
+    transversal area. Returns the internal-efforts at both extremal nodes.
+    '''
+    fi = np.zeros(2)
+    fi[0] = E * A * (Eu[0,0] + Eu[0,1])
+    fi[1] = E * A * (Eu[1,0] + Eu[1,1])
+    #Note that Eu is always diagonal. Its a truss, the kinematics is penalized.
+    return fi
+
+if __name__ == "__main__":
+
+    #
+    # Workflow
+    #
+    ###################################################################
+    ## 0) Definition of a mechanical state
+    ###################################################################
+    x , y , z = np.array([0,1]) , np.array([0,1]), np.array([0,1]) # This is the initial state
+    ###################################################################
+    ## 1) From the previous two-solutions, compute the new direction
+    ###################################################################
+    l = np.sqrt((x[0] - x[1])**2 + (y[0] - y[1])**2 + (z[0] - z[1])**2) # This is the element's length
+    cos_theta = get_cos(x,y,z,l)
+    print('the angles',type(cos_theta),cos_theta)
+    ###################################################################
+    ## 2) Evaluate the transformation matrix
+    ###################################################################
+    A = map_global_to_local(cos_theta[0],cos_theta[1],cos_theta[2])
+    print('transformation matrix')
+    print(A)
+    ###################################################################
+    ## 3) Convert the global force in 3d space to a 2d space.
+    ###################################################################
+    fg = 9.8 * np.array([0,0,-1,0,0,-1]) # This is a hand made external loading
+                                         # in our case gravity, negative z direction
+
+    fl = vec_global_to_local(A,fg)
+    print('global force',fg)
+    print('local force',fl)
+    ###################################################################
+    ## 4) Convert the global coordinates in 3d space to a 2d space.
+    ###################################################################
+    ug = np.array([x[0],x[1],y[0],y[1],z[0],z[1]])
+    ul = vec_global_to_local(A,ug)
+    print('global coords',ug)
+    print('local coords',ul)
+    ###################################################################
+    ## 5) Here starts the isoparametric evaluations. Integrate the external-forces
+    ###################################################################
+    npts = 1 # Number of Gauss points
+    int_fl = elemental_force(fl,ul,npts)
+
+    print('integrate force',int_fl)
+    ###################################################################
+    ## 5) Integrate the internal-forces
+    ###################################################################
+    E , A = 1 , 1
+    int_def = Green_Lagrange(npts,ul,2*ul)
+    print('The green lagrange strain tensor')
+    print(int_def)
+
+    int_force = elemental_int_force(int_def,E,A)
+
+    print('The internal forces')
+    print(int_force)
+
+    sys.exit()
+