assign values from data Array to GenericVector?

Question

Hi,
in the context of elastodynamics, i want to solve a problem of soil-structure interaction with some inclusions on the domain with different material properties, and i want to define these material properties as a DG Function using a MeshFunctions (type double), e.g., i tried this without succes:

using namespace dolfin;

  // Inclussion
  class Omega0 : public SubDomain
  {
    bool inside(const Array<double>& x, bool on_boundary) const
    {
      if (x[1] <= -0.75 && x[1] >= -1.75 && std::abs(x[0]) <= 1)
        return true;
      else 
        return false;
    }
  };

  // Mesh dimesions
  double b = 0.3, L = 2.5, h = 2.5, Lpml = 1;
  Point vertex_1(-(L + Lpml), 0);
  Point vertex_2(L + Lpml, -(h + Lpml));

  // Density of the mesh
  int nx = 30, ny = 0.5*nx;

int main()
{
  // Mesh
  RectangleMesh mesh(vertex_1, vertex_2, nx, ny, "left/right");

  // Create mesh functions over the cells
  MeshFunction<double> density(mesh, mesh.topology().dim());
  density = 1.0;     // Set all domain density equal to 1

  Omega0 Omega_0;
  Omega_0.mark(density, 2.0);    // Set density on Omega0 equal to 2

  DG0::FunctionSpace Vdg(mesh);
  Function rho(Vdg);

  // This is what i want to do (or something similar)
  rho.vector() = density.values();  // ¿GenericVector = Array?
}

it is possible to do something like this?

Answer 1

Hi, consider

  std::vector<double> values(rho.vector()->local_size(), 0);
  std::shared_ptr<const GenericDofMap> dofmap = Vdg.dofmap();
  for (CellIterator c(mesh); !c.end(); ++c)
  {
    values[(dofmap->cell_dofs(c->index()))[0]] = density[*c];
  }
  rho.vector()->set_local(values); 

Comments

Hi, thanks for your reply... there is any difference if i change

for (CellIterator c(mesh); !c.end(); ++c)
{
  values[(dofmap->cell_dofs(c->index()))[0]] = density[*c];
}

to

for (CellIterator c(mesh); !c.end(); ++c)
{
  values[c->index()] = density[*c];
}

?

---

Hi, the change you propose works fine here. However, it is based on an assumption that the cell order is the same as the dof order in the DG0 space. It is true now and it is unlikely that it will change in the future but you see that the dofmap makes the code more robust. Also it keeps you reminded that mesh entities and dofs of function spaces are in general enumerated in a different way, see e.g. here.

---

i got it... many thanks!