adding more plasticity equations (doesn't work yet)

elasticity_fufem.py

@@ -6,8 +6,9 @@ import dune.geometry
 import dune.grid
 import dune.functions
 
+from math import sqrt
 from dune.grid import DataType
-from dune.functions import defaultGlobalBasis, Power, Composite, Lagrange
+from dune.functions import defaultGlobalBasis, subspaceBasis, Power, Composite, Lagrange
 
 from fufemforms import *
 from utilities import *
@@ -18,6 +19,7 @@ def isNear(a,b):
 
 
 dimension = 2
+plasticDimension = int((dimension * (dimension + 1))/2 - 1) # This one here is the dimension of symmetric trace free matrices
 ansatzOrder = 2
 
 # Standard acceleration of gravity
@@ -52,40 +54,49 @@ def isDirichlet(x):
 def dirichletValues(x):
     return np.zeros(dimension)
 
-
+sqrtTwoHalf = sqrt(2.0)/2.0
 
 def globalAssembler(basis):
 
+    dispBasis = subspaceBasis(basis,0)
+    plasticBasis = subspaceBasis(basis,1)
     N = len(basis)
+    #print(N)
 
     # Mark all Dirichlet DOFs
     isConstrained = np.zeros( N, dtype=bool )
-    basis.interpolate(isConstrained,isDirichlet)
+    dispBasis.interpolate(isConstrained,isDirichlet)
 
     # Interpolate the boundary values
     constraintDOFValues = np.zeros( N )
-    basis.interpolate(constraintDOFValues,dirichletValues)
+    dispBasis.interpolate(constraintDOFValues,dirichletValues)
 
     # Identity matrix
     
     Id = dune.common.FieldMatrix_double_3_3([[1,0,0], [0,1,0], [0,0,1]]) if dimension==3 else dune.common.FieldMatrix_double_2_2([[1,0], [0,1]])
 
-
+    # This should be a range operator equation.
+    Epls = lambda h: sqrtTwoHalf* dune.common.FieldMatrix_double_2_2([[h[0], h[1]],[h[1],-h[0]]])
     # Symmetric gradient
     E = lambda w: symmetrize(grad(w))
 
     # Trial and test functions for variational problem
-    u = trialFunction(basis)
-    v = testFunction(basis)
+    u = trialFunction(basis, 0)
+    theta = trialFunction(basis, 1)
+    v = testFunction(basis, 0)
+    eta = testFunction(basis,1)
+
+    #help(theta)
 
     # Vector valued rhs coefficient function
     f = Coefficient(rhs, basis.gridView, rangeType=VectorType(dimension))
 
     # Use St. Venant-Kirchhoff material
-    sigma = 2*mu*E(u) + lambda_*(Id*trace(E(u)))
+    elasticStrains = E(u)-Epls(theta) 
+    sigma = 2*mu*(elasticStrains) + lambda_*(Id*trace(elasticStrains))
 
     # Bilinear form and rhs functional
-    a = integrate(dot(sigma, E(v)))
+    a = integrate(dot(sigma, E(v)-Epls(eta)))
     b = integrate(dot(f,v))
 
     # Assemble into matrix and vector
@@ -103,20 +114,38 @@ def globalAssembler(basis):
 # Create grid on a 100x1x1 beam
 grid = dune.grid.structuredGrid([0,0],[5,5],[40, 40])
 
-# Create a vector-valued nodal Lagrange FE basis
-basis = defaultGlobalBasis(grid, Power(Lagrange(order=ansatzOrder),exponent=dimension,blocked=False,layout="interleaved"))
+# Create a vector-valued nodal Lagrange FE basis for Displacements and plastic strains
+basis = defaultGlobalBasis(grid,Composite(Power(Lagrange(order=ansatzOrder),exponent=dimension,blocked=False,layout="interleaved"),Power(Lagrange(order=0),exponent=plasticDimension,blocked=False,layout="interleaved")))
 
-print("Dimension of FE space is "+str(len(basis)))
+numElements = grid.size(0)
+numDofs = len(basis)
+
+numElasticDofs = numDofs - plasticDimension * numElements
+numPlasticDofs = numDofs - numElasticDofs
+# Let's talk about the elastic problem
+displacementBasis = subspaceBasis(basis,0)
+plasticStrainBasis = subspaceBasis(basis,1)
 
+print(basis.size([1,4]))
+
+print("Number of Functionsspaces is "+str(dimension + plasticDimension))
+print("Number of Displacement Functionspaces are "+str(dimension))
+print("Number of plastic Strain Functionsspaces are "+str(plasticDimension))
+
+print("Dimension of FE space is "+str(len(basis)))
+print("Dimension of Displacement FE space should be "+str(numElasticDofs))
+print("Dimension of Plastic Strain FE space should be "+str(numPlasticDofs))
 # Compute A and b
-A,b = globalAssembler(basis)
+A,b = globalAssembler(basis=basis)
+
+#print(A.shape)
 
 # Solve linear system!
 x = scipy.sparse.linalg.spsolve(A, b)
 
-u = basis.asFunction(x)
+# u = basis.asFunction(x)
 
-vtk = grid.vtkWriter(0)
-u.addToVTKWriter("sol", vtk, DataType.PointVector)
-vtk.write("linear-elasticity-py")
+# vtk = grid.vtkWriter(0)
+#u.addToVTKWriter("sol", vtk, DataType.PointVector)
+#vtk.write("linear-elasticity-py")