This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

assign values from data Array to GenericVector?

0 votes

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?

asked Sep 28, 2015 by hernan_mella FEniCS Expert (19,460 points)
edited Oct 2, 2015 by hernan_mella

1 Answer

+1 vote
 
Best answer

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); 
answered Oct 3, 2015 by MiroK FEniCS Expert (80,920 points)
selected Oct 6, 2015 by hernan_mella

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!

...