Using trunk C++ API: cannot set MeshDomains

On 15 October 2011 21:51, Cian Wilson
<email address hidden> wrote:
> New question #174566 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/174566
>
> Hello,
>
> I'm trying (perhaps prematurely) to keep my code reasonably up to date with development trunk so was investigating the new MeshDomains class.

The more testing the better! You're also likely to get things fixed
sooner using the dev version.

> I'm reading a MeshFunction in from an xml file and trying to set the facet domains.
>
> Where previously I used:
>  dolfin::MeshFunction<dolfin::uint> edgeids(*mesh, filename);
>
>  boost::shared_ptr< dolfin::MeshFunction< dolfin::uint > > meshfuncedgeids =
>                           mesh->data().create_mesh_function("exterior_facet_domains");
>
>  *meshfuncedgeids = edgeids;
>
> I'm now trying:
>  dolfin::MeshFunction<dolfin::uint> edgeids(*mesh, filename);
>
>  const uint D = mesh->topology().dim();
>
>  mesh->domains().markers_shared_ptr(D-1).reset( new dolfin::MeshValueCollection<dolfin::uint>(edgeids) );
>
>  std::cerr << mesh->domains().markers(D-1).str(false) << std::endl;
>  if(!mesh->domains().facet_domains(*mesh))
>    dolfin::error("No facet_domains found.");
>
> which returns:
>  <MeshValueCollection of topological dimension 1 containing 0 values>
>  terminate called after throwing an instance of 'std::runtime_error'
>    what():  *** Error: No facet_domains found.
>
> Do you have any tips on how this should be done or is this still under development and likely to change?
>

It hopefully won't change. Can you send a complete (bit simple)
program, including the file that you're reading in.

Garth

> Many thanks for your help and for the great library!
> Cheers,
> Cian
>
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>

I can reproduce the problem in a much simpler case without actually reading in a mesh or mesh function (have also made it a bit more vebose):

main.cpp:
-------------------------------------------------------------
#include <dolfin.h>

class Side : public dolfin::SubDomain
{
    public:
    Side(const double &side) : side_(side)
    {
    }

    bool inside(const dolfin::Array<double>& x, bool on_boundary) const
    {
      return (std::fabs(x[0] - side_) < DOLFIN_EPS && on_boundary);
    }

  private:

    double side_;
};

int main()
{
  dolfin::UnitInterval mesh(10);

  Side left(0.0);
  Side right(1.0);

  dolfin::MeshFunction<dolfin::uint> edgeids(mesh, 0, 0);
  left.mark(edgeids, 1);
  right.mark(edgeids, 2);
  std::cout << "edgeids:" << std::endl;
  std::cout << edgeids.str(true) << std::endl;

  boost::shared_ptr< dolfin::MeshValueCollection<dolfin::uint> >
   edgecoll( new dolfin::MeshValueCollection<dolfin::uint>(edgeids) );
  std::cout << "edgecoll:" << std::endl;
  std::cout << edgecoll->str(false) << std::endl;

  std::cout << "makers(0) before assignment:" << std::endl;
  std::cout << mesh.domains().markers(0).str(false) << std::endl;

  mesh.domains().markers_shared_ptr(0) = edgecoll;

  std::cout << "makers(0) after assignment:" << std::endl;
  std::cout << mesh.domains().markers(0).str(false) << std::endl;
  if(!mesh.domains().facet_domains(mesh))
    dolfin::error("No facet_domains found.");
  else
    std::cout << "facet_domains found." << std::endl;

  return 0;
}
------------------------------------------------------------------

May I check if anyone can reproduce hitting the dolfin::error (the above main.cpp compiles with a CMakeLists.txt from a demo) or is this problem local to me? It's still affecting me when using 1.0-beta2.

Cheers,
Cian

I too hit the dolfin::error in dolfin-dev.

Your main.cpp file would work if

MeshDomains::markers_shared_ptr() (dolfin/mesh/MeshDomains.h)

    /// Get subdomain markers for given dimension (shared pointer version)
    boost::shared_ptr<MeshValueCollection<unsigned int> >
    markers_shared_ptr(uint dim);

returned a reference to a shared_ptr i.e.,

    /// Get subdomain markers for given dimension (shared pointer version)
    boost::shared_ptr<MeshValueCollection<unsigned int> >&
    markers_shared_ptr(uint dim);

I don't know if that's an acceptable solution?

Note that the MeshValueCollection you create from the MeshFunction 'edgeids' contains 20 values, while doing

  left.mark(*mesh.domains().markers_shared_ptr(0), 1, mesh);
  right.mark(*mesh.domains().markers_shared_ptr(0), 2, mesh);

or

  left.mark(mesh, 0, 1);
  right.mark(mesh, 0, 2);

results in mesh.domains().markers(0) having only 2 values. Don't know if there's a performance penalty involved using the former approach?

On Thursday October 27 2011 04:55:55 Kristian B. Ølgaard wrote:
> Question #174566 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/174566
>
> Kristian B. Ølgaard proposed the following answer:
> Your main.cpp file would work if
>
> MeshDomains::markers_shared_ptr() (dolfin/mesh/MeshDomains.h)
>
> /// Get subdomain markers for given dimension (shared pointer version)
> boost::shared_ptr<MeshValueCollection<unsigned int> >
> markers_shared_ptr(uint dim);
>
> returned a reference to a shared_ptr i.e.,
>
> /// Get subdomain markers for given dimension (shared pointer version)
> boost::shared_ptr<MeshValueCollection<unsigned int> >&
> markers_shared_ptr(uint dim);
>
> I don't know if that's an acceptable solution?

I think that defeats the benefit of a shared_ptr.

I see two other possible solutions.

  1) Call MeshDomains::init(dim), and grab MeshValueCollection for that dim
     and inititalize it with your MeshValueCollection (or MeshFunction
     directly)
  2) Using MeshDomain::set_mesh_value_collection()

Neither of the above mentioned methods work per today. I would vote that
functionality for 1 is added.

In general I think we need to make MeshValueCollection easier to work with.
Adding operator=(const {MeshFunction,MeshValueCollection}) is a start and that
would solve this problem.

Johan

> Note that the MeshValueCollection you create from the MeshFunction
> 'edgeids' contains 20 values, while doing
>
> left.mark(*mesh.domains().markers_shared_ptr(0), 1, mesh);
> right.mark(*mesh.domains().markers_shared_ptr(0), 2, mesh);
>
> or
>
> left.mark(mesh, 0, 1);
> right.mark(mesh, 0, 2);
>
> results in mesh.domains().markers(0) having only 2 values. Don't know if
> there's a performance penalty involved using the former approach?

Thanks Kristian and Johan.

I've only had a chance to play with Kritian's suggestions so far. The suggested reference method does indeed fix the problem but I understand if that's not the ideal solution. Will have a look at your suggested alternative modifications Johan.

I realise the example given was a bit contrived, as it doesn't make use of the sparse data storage of the MeshValueCollections but I was thinking of the case where dolfin-convert returns a MeshFunction that's read in from xml. I'm not aware of functionality that preserves the sparsity in this case - other than looping over the entries and assuming a 0 marker is equivalent to no id. Still, it's useful to know that the mark functionality is there for applying subdomains directly to the mesh domains.

Using mark directly on the mesh also fixes this particular example but I run into problems using that with a case that actually calls assemble (on ufc with interior facet integrals). In this case both the exterior_facet_domains and interior_facet_domains MeshFunctions get extracted from the facet domains at Assembler.cpp:131 and Assembler.cpp:139 respectively. Using the previous mesh data based method I would not have had an interior_facet_domain MeshFunction associated. This wouldn't be a problem except that the sparse data storage of the MeshValueCollection gets converted to dense MeshFunction data storage at MeshDomains.cpp:172 by initializing the MeshFunction to the max uint. This then falls foul of the check at Assembler.cpp:405 and none of my interior facets get assembled. Initializing the MeshFunction to 0 instead (at MeshDomains.cpp:172) fixes this but I can imagine that going wrong in some cases. Perhaps this is a separate bug report.

I've actually always been curious about the check at Assembler.cpp:405 (and its cell and exterior facet integral equivalents) because it means that if my mesh has domains attached, terms in my form that just use dx, ds or dS (rather than dx(id), etc.) don't get assembled, which isn't what I would have expected. Is this an intentional design decision? Just idle curiosity so perhaps this should be a separate question!

Cheers,
Cian

On 27 October 2011 18:20, Johan Hake
<email address hidden> wrote:
> Question #174566 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/174566
>
> Johan Hake proposed the following answer:
> On Thursday October 27 2011 04:55:55 Kristian B. Ølgaard wrote:
>> Question #174566 on DOLFIN changed:
>> https://answers.launchpad.net/dolfin/+question/174566
>>
>> Kristian B. Ølgaard proposed the following answer:
>> Your main.cpp file would work if
>>
>> MeshDomains::markers_shared_ptr() (dolfin/mesh/MeshDomains.h)
>>
>>     /// Get subdomain markers for given dimension (shared pointer version)
>>     boost::shared_ptr<MeshValueCollection<unsigned int> >
>>     markers_shared_ptr(uint dim);
>>
>> returned a reference to a shared_ptr i.e.,
>>
>>     /// Get subdomain markers for given dimension (shared pointer version)
>>     boost::shared_ptr<MeshValueCollection<unsigned int> >&
>>     markers_shared_ptr(uint dim);
>>
>> I don't know if that's an acceptable solution?
>
> I think that defeats the benefit of a shared_ptr.

Yes, me too.

> I see two other possible solutions.
>
>  1) Call MeshDomains::init(dim), and grab MeshValueCollection for that dim
>     and inititalize it with your MeshValueCollection (or MeshFunction
>     directly)
>  2) Using MeshDomain::set_mesh_value_collection()

Or simply let MeshDomain::_markers be public, to me this usage seems
relatively low level and only advanced users should need it, but I
could be very wrong of course.

> Neither of the above mentioned methods work per today. I would vote that
> functionality for 1 is added.
>
> In general I think we need to make MeshValueCollection easier to work with.
> Adding operator=(const {MeshFunction,MeshValueCollection}) is a start and that
> would solve this problem.

If this was added, I think we could get away with doing

*mesh.domains().markers_shared_ptr(0) = edgecoll;

in the code.

Kristian

> Johan
>
>> Note that the MeshValueCollection you create from the MeshFunction
>> 'edgeids' contains 20 values, while doing
>>
>>   left.mark(*mesh.domains().markers_shared_ptr(0), 1, mesh);
>>   right.mark(*mesh.domains().markers_shared_ptr(0), 2, mesh);
>>
>> or
>>
>>   left.mark(mesh, 0, 1);
>>   right.mark(mesh, 0, 2);
>>
>> results in mesh.domains().markers(0) having only 2 values. Don't know if
>> there's a performance penalty involved using the former approach?
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>

I'd suggest reporting a bug with a very simple piece of complete code that demonstrates the issue. That way we won't lose track of it.

On Thursday October 27 2011 12:15:51 Kristian B. Ølgaard wrote:
> Question #174566 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/174566
>
> Kristian B. Ølgaard proposed the following answer:
> On 27 October 2011 18:20, Johan Hake
>
> <email address hidden> wrote:
> > Question #174566 on DOLFIN changed:
> > https://answers.launchpad.net/dolfin/+question/174566
> >
> > Johan Hake proposed the following answer:
> >
> > On Thursday October 27 2011 04:55:55 Kristian B. Ølgaard wrote:
> >> Question #174566 on DOLFIN changed:
> >> https://answers.launchpad.net/dolfin/+question/174566
> >>
> >> Kristian B. Ølgaard proposed the following answer:
> >> Your main.cpp file would work if
> >>
> >> MeshDomains::markers_shared_ptr() (dolfin/mesh/MeshDomains.h)
> >>
> >> /// Get subdomain markers for given dimension (shared pointer
> >> version) boost::shared_ptr<MeshValueCollection<unsigned int> >
> >> markers_shared_ptr(uint dim);
> >>
> >> returned a reference to a shared_ptr i.e.,
> >>
> >> /// Get subdomain markers for given dimension (shared pointer
> >> version) boost::shared_ptr<MeshValueCollection<unsigned int> >&
> >> markers_shared_ptr(uint dim);
> >>
> >> I don't know if that's an acceptable solution?
> >
> > I think that defeats the benefit of a shared_ptr.
>
> Yes, me too.
>
> > I see two other possible solutions.
> >
> > 1) Call MeshDomains::init(dim), and grab MeshValueCollection for that
> > dim and inititalize it with your MeshValueCollection (or MeshFunction
> > directly)
> > 2) Using MeshDomain::set_mesh_value_collection()
>
> Or simply let MeshDomain::_markers be public, to me this usage seems
> relatively low level and only advanced users should need it, but I
> could be very wrong of course.

Sure, that is also a way. Is they private for a reason?

> > Neither of the above mentioned methods work per today. I would vote that
> > functionality for 1 is added.
> >
> > In general I think we need to make MeshValueCollection easier to work
> > with. Adding operator=(const {MeshFunction,MeshValueCollection}) is a
> > start and that would solve this problem.
>
> If this was added, I think we could get away with doing
>
> *mesh.domains().markers_shared_ptr(0) = edgecoll;

Exactly!

Johan

> in the code.
>
> Kristian
>
> > Johan
> >
> >> Note that the MeshValueCollection you create from the MeshFunction
> >> 'edgeids' contains 20 values, while doing
> >>
> >> left.mark(*mesh.domains().markers_shared_ptr(0), 1, mesh);
> >> right.mark(*mesh.domains().markers_shared_ptr(0), 2, mesh);
> >>
> >> or
> >>
> >> left.mark(mesh, 0, 1);
> >> right.mark(mesh, 0, 2);
> >>
> >> results in mesh.domains().markers(0) having only 2 values. Don't know if
> >> there's a performance penalty involved using the former approach?
> >
> > --
> > You received this question notification because you are a member of
> > DOLFIN Team, which is an answer contact for DOLFIN.

On 10/27/2011 06:55 PM, Johan Hake wrote:
> Your question #174566 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/174566
>
> Johan Hake proposed the following answer:
> On Thursday October 27 2011 12:15:51 Kristian B. Ølgaard wrote:
>> Question #174566 on DOLFIN changed:
>> https://answers.launchpad.net/dolfin/+question/174566
>>
>> Kristian B. Ølgaard proposed the following answer:
>> On 27 October 2011 18:20, Johan Hake
>>
>> <email address hidden> wrote:
>>> Question #174566 on DOLFIN changed:
>>> https://answers.launchpad.net/dolfin/+question/174566
>>>
>>> Johan Hake proposed the following answer:
>>>
>>> On Thursday October 27 2011 04:55:55 Kristian B. Ølgaard wrote:
>>>> Question #174566 on DOLFIN changed:
>>>> https://answers.launchpad.net/dolfin/+question/174566
>>>>
>>>> Kristian B. Ølgaard proposed the following answer:
>>>> Your main.cpp file would work if
>>>>
>>>> MeshDomains::markers_shared_ptr() (dolfin/mesh/MeshDomains.h)
>>>>
>>>> /// Get subdomain markers for given dimension (shared pointer
>>>> version) boost::shared_ptr<MeshValueCollection<unsigned int> >
>>>> markers_shared_ptr(uint dim);
>>>>
>>>> returned a reference to a shared_ptr i.e.,
>>>>
>>>> /// Get subdomain markers for given dimension (shared pointer
>>>> version) boost::shared_ptr<MeshValueCollection<unsigned int> >&
>>>> markers_shared_ptr(uint dim);
>>>>
>>>> I don't know if that's an acceptable solution?
>>> I think that defeats the benefit of a shared_ptr.
>> Yes, me too.
>>
>>> I see two other possible solutions.
>>>
>>> 1) Call MeshDomains::init(dim), and grab MeshValueCollection for that
>>> dim and inititalize it with your MeshValueCollection (or MeshFunction
>>> directly)
>>> 2) Using MeshDomain::set_mesh_value_collection()
>> Or simply let MeshDomain::_markers be public, to me this usage seems
>> relatively low level and only advanced users should need it, but I
>> could be very wrong of course.
> Sure, that is also a way. Is they private for a reason?
>
>>> Neither of the above mentioned methods work per today. I would vote that
>>> functionality for 1 is added.
>>>
>>> In general I think we need to make MeshValueCollection easier to work
>>> with. Adding operator=(const {MeshFunction,MeshValueCollection}) is a
>>> start and that would solve this problem.
>> If this was added, I think we could get away with doing
>>
>> *mesh.domains().markers_shared_ptr(0) = edgecoll;
> Exactly!
>
> Johan
>

I've implemented an operator=(const MeshFunction) and it works. In fact
it's even simpler as you don't have to use the shared_ptr interface
anymore, just:
mesh.domains().markers(0) = edgeids;

Cheers,
Cian

>> in the code.
>>
>> Kristian
>>
>>> Johan
>>>
>>>> Note that the MeshValueCollection you create from the MeshFunction
>>>> 'edgeids' contains 20 values, while doing
>>>>
>>>> left.mark(*mesh.domains().markers_shared_ptr(0), 1, mesh);
>>>> right.mark(*mesh.domains().markers_shared_ptr(0), 2, mesh);
>>>>
>>>> or
>>>>
>>>> left.mark(mesh, 0, 1);
>>>> right.mark(mesh, 0, 2);
>>>>
>>>> results in mesh.domains().markers(0) having only 2 values. Don't know if
>>>> there's a performance penalty involved using the former approach?
>>> --
>>> You received this question notification because you are a member of
>>> DOLFIN Team, which is an answer contact for DOLFIN.

On 10/27/2011 06:15 PM, Garth Wells wrote:
> Your question #174566 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/174566
>
> Garth Wells posted a new comment:
> I'd suggest reporting a bug with a very simple piece of complete code
> that demonstrates the issue. That way we won't lose track of it.
>
I've spun off two bug reports - one for the original topic and another
for the failure to integrate internal facets.

See https://code.launchpad.net/~ldeo-magma/dolfin/devfixes/+merge/80709

On Thu, Oct 27, 2011 at 05:31:05PM -0000, Cian Wilson wrote:
> Question #174566 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/174566
>
> Cian Wilson posted a new comment:
> Thanks Kristian and Johan.
>
> I've only had a chance to play with Kritian's suggestions so far. The
> suggested reference method does indeed fix the problem but I understand
> if that's not the ideal solution. Will have a look at your suggested
> alternative modifications Johan.
>
> I realise the example given was a bit contrived, as it doesn't make use
> of the sparse data storage of the MeshValueCollections but I was
> thinking of the case where dolfin-convert returns a MeshFunction that's
> read in from xml. I'm not aware of functionality that preserves the
> sparsity in this case - other than looping over the entries and assuming
> a 0 marker is equivalent to no id. Still, it's useful to know that the
> mark functionality is there for applying subdomains directly to the mesh
> domains.
>
> Using mark directly on the mesh also fixes this particular example but I
> run into problems using that with a case that actually calls assemble
> (on ufc with interior facet integrals). In this case both the
> exterior_facet_domains and interior_facet_domains MeshFunctions get
> extracted from the facet domains at Assembler.cpp:131 and
> Assembler.cpp:139 respectively. Using the previous mesh data based
> method I would not have had an interior_facet_domain MeshFunction
> associated. This wouldn't be a problem except that the sparse data
> storage of the MeshValueCollection gets converted to dense MeshFunction
> data storage at MeshDomains.cpp:172 by initializing the MeshFunction to
> the max uint. This then falls foul of the check at Assembler.cpp:405
> and none of my interior facets get assembled. Initializing the
> MeshFunction to 0 instead (at MeshDomains.cpp:172) fixes this but I can
> imagine that going wrong in some cases. Perhaps this is a separate bug
> report.
>
> I've actually always been curious about the check at Assembler.cpp:405
> (and its cell and exterior facet integral equivalents) because it means
> that if my mesh has domains attached, terms in my form that just use dx,
> ds or dS (rather than dx(id), etc.) don't get assembled, which isn't
> what I would have expected. Is this an intentional design decision?
> Just idle curiosity so perhaps this should be a separate question!

I don't know if this question was answered, but if domains are given
as input to the assembler in any way (there are several ways for it to
enter), then DOLFIN assumes that data provides a numbering of the
subdomain so then that is used. And dx is short for dx(0) in the form
language. At some point, we discussed letting dx mean "all subdomains"
but it would require some redesign of the UFC interface.

--
Anders

Thanks Anders, it was mostly curiosity as dx meaning dx(0) sometimes trips me up.

The question did spin off a bug report however about passing the domains to the assembler using MeshDomain data:
https://bugs.launchpad.net/dolfin/+bug/882910
The assumption that dS means dS(0) and the setting of the the MeshFunction to max uint when converting from a MeshValueCollection for the interior (as well as for the outer) facet domains means that the behaviour changed between 1.0-beta (MeshFunctions beneath MeshData) and 1.0-beta2 (MeshValueCollections beneath MeshDomains). Essentially interior facets are no longer assembled in the given example. I can still reproduce this with dev though some of the checks at the end of main.cpp need modifying to take account of vector() now returning a boost shared_ptr.