dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Wed, Aug 13, 2008 at 01:18:32PM -0500, Catherine Micek wrote:
> Hi,
> 
> I'm a (relatively) new user to the Fenics system and have a  
> question.  I'm not sure if this is the correct place to send this,  
> but I couldn't find anything
> in the repository or the manuals (perhaps there's a relevant example  
> somewhere I missed ... ?), so I ended up here.
> 
> I'm using Fenics to study the equations of linear elasticity (in a  
> mixed Stokes formulation), which allows me to directly compute a  
> displacement U and
> pressure P.  In addition to these two variables, I am interested in  
> the stress of the system. I have been able to obtain the stress by  
> defining a suitable
> projection:
> 
> 	stresselement = VectorElement("Discontinuous Lagrange", "triangle",  
> 1, 9)
> 	S = TrialFunction(stresselement)
>     	T = TestFunction(stresselement)
> 
> 	a = dot(S, T)*dx
> 	# Sigma is the symmetric stress computed from the displacement, U
>     	L = dot(Sigma, t)*dx
> 
> and then solving the resulting linear system for S.  Next, I would  
> like to do some calculations with the stresses; for example,  
> integrate the stress along a boundary.
> To do this, I need the stress vector, s = Sn, where n is the normal  
> on the boundary.  However, I am unclear how to do this, as I have  
> defined the stress as a 9-component
> vector instead of a 3x3 matrix (and n will be defined as a 3  
> component vector if I use the FacetNormal function).
> 
> Any suggestions on how to handle this?  Thanks!
> 
> Katy

We frequently compute normal and shear stresses on the boundary for
Navier-Stokes from a computed velocity field U and pressure P.

I've attached the code we use. See if you can make sense of it.

-- 
Anders

# Computation of normal and shear stresses on
# the boundary for a given solutioin (u, p)
#
# Anders Logg and Kristian Valen-Sendstad 2008-06-12 -- 2008-06-21

from dolfin import *

def epsilon(u):
    "Return strain-rate tensor"
    return 0.5*(grad(u) + transp(grad(u)))

class Stress:
    "Computation of stress for given u, p"

    def __init__(self, u, p, nu, mesh):
        "Initialize stress"

        # Choose shape
        if mesh.topology().dim() == 2:
            shape = "triangle"
        else:
            shape = "tetrahedron"

        # Compute stress tensor
        sigma = mult(2.0*nu, epsilon(u)) - mult(p, Identity(len(u)))

        # Compute stress on surface
        n = FacetNormal(shape, mesh)
        F = mult(sigma, -n)

        # Compute normal and tangential components
        Fn = dot(F, n)       # scalar-valued
        Ft = F - mult(Fn, n) # vector-valued

        # Integrate against piecewise constants on the boundary
        scalar  = FiniteElement("DG", shape, 0)
        vector  = VectorElement("DG", shape, 0)
        scaling = InvFacetArea(shape, mesh)
        v = TestFunction(scalar)
        w = TestFunction(vector)
        self.Ln = scaling*v*Fn*ds
        self.Lt = scaling*dot(w, Ft)*ds

        # Create functions
        self.Fn = Function(scalar, mesh, Vector())
        self.Ft = Function(vector, mesh, Vector())

        self.mesh = mesh

    def __call__(self):
        "Return normal and shear stress (Fn, Ft)"

        print "Computing stress"
        
        # Assemble vectors
        assemble(self.Ln, self.mesh, tensor=self.Fn.vector())
        assemble(self.Lt, self.mesh, tensor=self.Ft.vector())

        # FIXME: Cannot extract boundary function yet, because of problems
        # FIXME: with embedding element in higher dimension.
        #cell_map = self.boundary.data().meshFunction("cell map")
        #for boundary_cell in cells(self.boundary):
        #    mesh_facet = Facet(self.mesh, cell_map(boundary_cell))
        #    print mesh_facet.index()
        #    for mesh_cell in cells(mesh_facet):
        #        print bn[mesh_cell.index()]
        #        print boundary_cell.index(), xn.size()

        return (self.Fn, self.Ft)

Attachment: signature.asc
Description: Digital signature