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

Steffen Müthing steffen.muething at iwr.uni-heidelberg.de
Wed Jun 13 11:00:18 CEST 2018 


> Am 13.06.2018 um 10:51 schrieb Michael Wenske <m_wens01 at wwu.de>:
> 
> 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.

Good catch! That’s very much a bug, yes. If you want to file it, just register
at gitlab.dune-project.org and go to https://gitlab.dune-project.org/pdelab/dune-pdelab.

Cheers
Steffen

> 
> 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
>> 
> 
> 
> _______________________________________________
> dune-pdelab mailing list
> dune-pdelab at lists.dune-project.org
> https://lists.dune-project.org/mailman/listinfo/dune-pdelab

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 488 bytes
Desc: Message signed with OpenPGP
URL: <https://lists.dune-project.org/pipermail/dune-pdelab/attachments/20180613/882153ae/attachment.sig>

More information about the dune-pdelab mailing list