[dune-pdelab] Dune PDELab (2.5) Constraints questions

[Michael Wenske]

Hi!

I have a hunch that the problem lies within the struct VolumeConstraints
in constraints.hh. Within this struct there are these two lines:

        // allocate local constraints map
        CL cl;

which seem to have the effect that when the constraints assembler
calls the applyToTree method on the volume constraints, the transformation
which is handed down to the user-implemented VolumeConstraints is
disregarded
and overwritten with an empty one. The above line overwrites the trafo
handed into the constructor, if I am correct.

The other constraints structs SkeletonConstraints, ProcessorConstraints,
BoundaryConstraints etc
do not have a comparable statement, and if I comment it out everything
works as intended and the
user is free to implement his/her volume constraints as needed.

I also see that same line in the dune-pdelab/2.6 branch within git. If
this is correct,
what is the proper procedure to file a bug report / merge request?

thanks,

Michael

On 12.06.2018 16:17, Christian Engwer wrote:
> Hi Michael.
>
> I have somewhere a separate assembler to constrain certain parts of
> the volume (we used this for a quick and dirty implementation of
> surface FEM...). You can also ask Sebastian, he should have the
> implementation.
>
> Best
> Christian
>
> On Tue, Jun 12, 2018 at 01:18:41PM +0200, Steffen Müthing wrote:
>> Hi Michael,
>>
>>> Am 12.06.2018 um 11:07 schrieb Michael Wenske <m_wens01 at wwu.de>:
>>>
>>> Dear all!
>>>
>>> I am struggling with the Constraints implementation and the assembly of
>>> such in Dune::PDELab.
>>> I want to constrain all degrees of freedom within a subvolume(tumor) in
>>> the inside of my domain(brain).
>>>
>>> My degrees of freedom are an the edges of a quadrilateral Grid and are
>>> assembled, via the alpha_volume(...)
>>> methods in the local operator. I've studied files like
>>> pdelab/constraints/conforming.hh and their interface seems to imply that
>>> only fluxes over an -interface-  are
>>> to be dirichlet constrained. Now theoretically there should be no
>>> difference wether a dof is restricted by a loop
>>> which is iterating over the interfaces, or if it is iterating the volumes.
>> The currently implemented constraints assemblers all work on intersections
>> because that is what you need most of the time. No one has needed a volume
>> constraints assembler until now, so there isn’t one in PDELab.
>>
>>> Question 1: can I use the skeleton constraints and volume constraints
>>> interchangeably as long as i
>>> hit the right dof's and set their transformation to empty?
>> But you are right, you can write your own constraints assembler (like you did)
>> and have that constrain DOFs on a cell instead of DOFs on an intersection.
>>
>>> I see that there is an interface to pass an object to the constraints
>>> assembler which might be used
>>> to let the user define the constraints via a simple lambda. something like:
>>>
>>> using CON = Dune::PDELab::OverlappingConformingDirichletConstraints;
>>> using CC = typename GFS::template
>>> ConstraintsContainer<RangeFieldType>::Type;
>>> auto bc_lambda = [&](const auto& i,const auto& x) {return
>>> modelLOP_p->b(i,x);};
>>> auto b = Dune::PDELab::makeBoundaryConditionFromCallable(*gv_p,bc_lambda);
>>>
>>> Dune::PDELab::constraints(b,*gfs_p, *cc_p);
>>>
>>> My problem is that the above function will derive from
>>> DirichletConstraintsParameters, again
>>> being restricted to the outer boundary and not the inside. On a wider
>>> scale I see the ConstraintsParameters
>>> interfaces which I take are there to define the type of the object and
>>> signature of the methods contained in b.
>>> This object is then passed all the way down to the actual method
>>> deciding isDirichlet or isNeumann etc.
>>>
>>> 2. Is the XXXParameters concept to have the user write these
>>> interfaces...and then implement the constraints class using it him/herself?
>>> How would I use the Parameter Interface so that my constraints would
>>> work with existing implementations?
>> The existing parameter classes are written for the more common concept of boundary
>> constraints, so they will not work in your case. So it’s like you guessed: You have to write
>> your own constraints assembler, and you can design the parameter interface any way
>> you want.
>>
>>
>>> I did'nt give up yet. I wrote my own constraints class. It's a little
>>> class fulfilling (afaik) the interface for the
>>> VolumeConstraints given in pdelab/constraints/common/constraints.hh (see
>>> file)
>>>
>>> it will get an elementGeometry from the assembler, decide if this dof is
>>> to be restricted or not by the position of
>>> the dof's (isTumor(pos)) and then set the transformation to empty. I can
>>> debug right into the line where the
>>> transformation is set to empty.
>>>
>>>             // check if this corner-dof is supposed to be constrained
>>>             if (medData_p->isTumor(pos)) trafo[lfs.dofIndex(i)] = empty;
>>>
>>> I can not grasp where i went wrong, and I guess there is some conceptual
>>> thing
>>> i am missing about the way these interfaces are written. When i run the
>>> program it tells me that there are -zero-
>>> constrained dof's in my gridfunctionspace even though, when I debug the
>>> assembly procedure seems to take
>>> place correctly. Also the way I now do it, I do not use the
>>> ConstraintsParameter Object passed into my class at all.
>> That’s strange, looking at the attached code, that should do the right thing. What exactly
>> are you doing in your main driver function? Are you sure that the isTumor() call actually returns
>> true?
>>
>>> 3. If anybody could give a quick overview about the -idea- behind the
>>> setup of constraints in pdelab , or what i am currently missing, i would
>>> be most
>>> grateful. What would be the cleanest way for a user to implement the
>>> constriction of dof's at different positions in the domain?
>> I think you already managed to figure it out pretty much on your own… :-)
>>
>>> I would also be gratefull for hits to any code examples where people
>>> have implemented interior boundary conditions before.
>> As I said: I’m personally not aware of such examples.
>>
>> Cheers
>> Steffen
>>
>>> thank you in advance,
>>>
>>> Michael
>>>
>>> <SegmentationConstraints.hh>_______________________________________________
>>> dune-pdelab mailing list
>>> dune-pdelab at lists.dune-project.org
>>> https://lists.dune-project.org/mailman/listinfo/dune-pdelab
>
>
>> _______________________________________________
>> dune-pdelab mailing list
>> dune-pdelab at lists.dune-project.org
>> https://lists.dune-project.org/mailman/listinfo/dune-pdelab
>