How to compute gradient of displacement and other functions in ex2.cpp

[GoogleCodeExporter]

I am new to MFEM. In example ex2.cpp I would like to compute gradient of 
displacement field u and then stress sigma(u), strain eps(u), strain energy U 
and elementwise strain energy Ue. The functional form of these quantities as a 
function of u is given by

eps(u) = (grad(u) + grad(u)^T)/2

sigma(u) = 2*mu*eps + lambda*trace(eps) I

U = lmbda/2.0*(trace(eps))^2 + mu*trace(eps^2) 

Ue= \int_(Omega)(U/vol *v dOmega) , where vol is cell volume and v are trial 
functions 

Can you provide an example or illustration for this?
Also, how is it possible to apply a point load at specific location on free end 
of the cantilever beam (ex2.cpp)?

Original issue reported on code.google.com by yogeshpa...@gmail.com on 17 Mar 2014 at 10:36

[GoogleCodeExporter]

One way to compute fields derived from a solution is to define a special
sub-class of the Coefficient class and then project it on a given GridFunction.
Here is an example:

class MyCoefficient : public Coefficient
{
private:
   GridFunction &u;
   Coefficient &lambda, μ
   DenseMatrix eps, sigma;

public:
   MyCoefficient(GridFunction &_u, Coefficient &_lambda, Coefficient &_mu)
      : u(_u), lambda(_lambda), mu(_mu) { }
   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
   {
      u.GetVectorGradient(T, eps);  // eps = grad(u)
      eps.Symmetrize();             // eps = (1/2)*(grad(u) + grad(u)^t)
      double l = lambda.Eval(T, ip);
      double m = mu.Eval(T, ip);
      sigma.Diag(l*eps.Trace(), eps.Size()); // sigma = lambda*trace(eps)*I
      sigma.Add(2*m, eps);          // sigma += 2*mu*eps

      return sigma(0,0); // return sigma_xx
   }
   virtual void Read(istream &in) { }
   virtual ~MyCoefficient() { }
};

Then inside the main function in ex2.cpp we can project this coefficient to a
GridFunction and save it to a file:

   {
      // A. Define a finite element space for post-processing the solution. We
      //    use a discontinuous space of the same order as the solution.
      L2_FECollection pp_fec(p, dim);
      FiniteElementSpace pp_fespace(mesh, &pp_fec);
      GridFunction pp_field(&pp_fespace);

      // B. Project the post-processing coefficient defined above to the
      //    'pp_field' GridFunction.
      MyCoefficient pp_coeff(x, lambda_func, mu_func);
      pp_field.ProjectCoefficient(pp_coeff);

      // C. Save the 'pp_field'.
      ofstream pp_ofs("postproc.gf");
      pp_ofs.precision(8);
      pp_field.Save(pp_ofs);
   }


Element-wise strain energies can be computed using a similar approach: in this
case the coefficient can be used to define a LinearForm (on a piece-wise
constant FE space) with one DomainLFIntegrator using a strain energy density
coefficient. After assembling, the LinearForm will give you the element-wise
energies.


Regarding a point load, one way to achieve this is by using a sub-class of
VectorCoefficient, something like this:

class PointLoadCoefficient : public VectorCoefficient
{
private:
   Vector center, force, point;

public:
   PointLoadCoefficient(const Vector &_center, const Vector &_force)
      : VectorCoefficient(_center.Size()), center(_center), force(_force) { }
   virtual void Eval(Vector &V, ElementTransformation &T,
                     const IntegrationPoint &ip)
   {
      T.Transform(ip, point);
      V.SetSize(vdim);
      if (point.DistanceTo(center) < 1e-12)
         V = force;
      else
         V = 0.0; // set all components to zero
   }
   virtual ~PointLoadCoefficient() { }
};

Projecting this coefficient on a GridFunction (defined on 'fespace' in ex2.cpp)
and adding it to the load vector should do the job.

Original comment by veselin@gmail.com on 22 Mar 2014 at 8:58

[GoogleCodeExporter]

Thank you!

Original comment by yogeshpa...@gmail.com on 23 Apr 2014 at 10:09

[GoogleCodeExporter]

Original comment by tzanio on 8 Dec 2014 at 10:00

• Changed state: Fixed

[fsabet3]

Thank you for the response. Why a discontinuous space (L2_FECollection) is used, while in example 2 a continuous space (H1_FECollection) was used?

[tzanio]

The gradient of H1 field is in L2.