How to assign different material properties to subdomains in C++ for a 3D/2D Gmsh?

[Chaitanya Raj Goyal]

Hello,

Suppose I create a 3D/ 2D mesh in gmsh. For Eg. Let's stick to a 3D mesh - a cuboid divided in 3 layers of different thickness like a cake. Each layer has 3 material properties - Young's Modulus, Poisson's ratio and Mass Density.

How can one define the individual material properties for each subdomain and use to solve a problem. I am looking at a non linear elastodynamic problem given here.

I understand how it can be done in python for a mesh created in FEniCS and also in Gmsh. 

Can someone please post a sample code or explain how to do it in C++ using Gmsh and if possible, also a fenics inbuilt mesh. I can post a sample Gmsh .geo file if that is required.

related to an answer for: How to set diferent materials on difrent parts of a mesh?

Thank you very much,

[David Bernstein]

Hi Chaitanya, This is definitely doable, though a bit complicated. In Gmsh you need to create physical volumes that correspond to each layer. The volume numbers are reflected by values in the .msh file that Gmsh creates. These are the material numbers for each cell. You then write a C++ function that reads the .msh file, extracts those numbers and creates a dolfin::MeshFunction<size_t>. The Young's modulus is then a dolfin::Function that uses this MeshFunction and a table of values for the Young's modulus for each material in its eval function.

I have done some FEniCS elasticity simulations using Gmsh in exactly this way and it works fine. Does require a fair amount of coding though.

Dave

[Chaitanya Raj Goyal]

Dear Dr. Bernstein,

Thank you very much for your email. I have successfully used the Gmsh physical tags to solve a dynamic elasticity problem with different material properties in python with FEniCS. I can also do it for a FEniCS inbuilt mesh in python. I think the problem is lack of coding experience in C++.

Right now I am experimenting with the following: (It works perfectly when run in serial, but fails in parallel)

 

auto V0 = std::make_shared<DG0::FunctionSpace>(mesh);
 
  MeshFunction<std::size_t> subdomains(mesh, mesh.topology().dim(), true);


  Omega0 Omega_0;            // This subdomain is defined as (x[1] <= 0.5) for a Gmsh square (1 x 1)
  Omega_0.mark(subdomains, 0);


  Omega1 Omega_1;              // This subdomain is defined as (x[1] >= 0.5) for a Gmsh square (1 x 1)
  Omega_1.mark(subdomains, 1);


  // Material parameters -->  mass density
  Function rho(V0);
  Array<double> rho_values(2);

  rho_values[0] = 0.5;
  rho_values[1] = 2;


  std:size_t cell_no, subdomain_no;
  for (CellIterator c(mesh); !c.end(); ++c)
  {
   cell_no = c->index();                                 // Cell number
   subdomain_no = subdomains.values()[cell_no];          // Number subdomain identifier that owns the cell number " cell_no "
   rho.vector()->setitem(cell_no, rho_values[subdomain_no]);  


   std::cout << "" << std::endl;
   std::cout << "- sub_domain_number: " << subdomains.values()[cell_no] << std::endl;
   std::cout << "- cell_number: " << cell_no << std::endl;
   std::cout << "- rho_in_cell: " << rho.vector()->getitem(cell_no) << std::endl;
  }


  plot(mesh);
  interactive();


  plot(subdomains);
  interactive();  


  plot(rho);
  interactive();  

Another issue with this approach is that I am not using Gmsh tags and this kind of subdomain division will become difficult for complicated geometry. If your elasticity code was in C++ and is not copyrighted, please consider sharing it with associated geo/msh file. That will be very useful.

Thanks a lot,

Sincerely,

Chaitanya

[David Bernstein]

I think the problem with your code in parallel may be the loop where you iterate over cells. This will not work in parallel since each process has only a subset of the cells.

Unfortunately our code is closed source so I can't share it. I would be happy to answer other questions though.

Dave

On Thursday, August 18, 2016 at 11:43:25 AM UTC-7, Chaitanya Raj Goyal wrote:

[Chaitanya Raj Goyal]

Dear All,

So, I am not pursuing the Gmsh tags with multiple materials issue in C++ anymore. I am just looking at a simple inbuilt fenics mesh or Gmsh, and trying to implement 'expression' with different values using if-else statement to define properties in different regions of the mesh. I tested this approach successfully with a simple C++ program where I just plotted the changing parameter (property) with the mesh. It worked in serial and parallel. Then I tried a Poisson problem. It again worked in serial and parallel. I have attached both these folders if you feel interested to run them and see for yourself.

Now, the problem is that I am trying to use Fenics Solid Mechanics app (link). In this library, you have the option to call 2 physical models in your code. A vonMises model or a DruckerPrager model. Currently, when I use vonMises model, and I try to compile the Makefile to create a demo, I get an error. I have attached the file - error.png. You can observe, the error is related to the vonMises class. So, in the FSM library, there are 2 files - VonMises.cpp and VonMises.h. If you look at these files, in the first line you' ll see that the vonMises class intakes the parameter 'E' (one of the parameters I divided) as a 'double'. 

class VonMises : public PlasticityModel { public: VonMises(double E, double nu, double yield_stress, double hardening_parameter);

However, not I am not defining it as a 'double' anymore in my code. Earlier, I used to define it in my code as a double:

double E = 1500.0;

and now I am defining it as,

class Young : public Expression
{
  void eval(Array<double>& values, const Array<double>& x) const
  {
    if (x[1] >= 0.5 )
    {
      values[0] = 27000.0;
    }
    else
    {
      values[0] = 15000.0;
    }
  }
};
Young E_v;

auto V = std::make_shared<e::FunctionSpace>(mesh);

auto E = std::make_shared<Function>(V);

E->interpolate(E_v);

The vonMises.h file also calls PlasticityModel.h in the beginning. Even that file identifies 'E' as a double and then uses it to compute const double 'mu' and const double 'lambda' for the elastic tangent. Can anyone please recommend what can be done about this? Should I remove double from the library and install the library again or change something in the code itself, so that it will be accepted by the FSM library. 

Please help with this issue.

Thanks a lot,

Chaitanya