How to store Constraints

[philbucher]

I am making this issue to continue the discussion started in #841

In a nutshell, #841 adds a BuilderAndSolver with Constraints and now the main discussion is where and how to store the constraints

I am making this issue bcs the discussion started in #841 but I think that this is interesting for others as well

[adityaghantasala]

@pooyan-dadvand @RiccardoRossi Following your last two comments in #841 :

I completely see your point in having a stable and long lasting solution for constraints in Kratos. Now given that we should also take care that this implementation/interface is independent of applications/physics.

Now coming to your comment saying

@adityaghantasala the point we are rising with @pooyan-dadvand is that a nice solution would be to treat the ConstraintEquation the same way as elements and conditions, so that a constraint equation can be added to a submodelpart and it will be automatically present in the parent modelparts. This would make the behaviour consistent with elements and conditions and would be a clean replacement for the ConstraintContainer (i don't recall the name now)

doing this would imply adding AddConstraint/AddConstraints as well as RemoveConstraints interface, as well as a few other methods that exist for Elements and Conditions.

This will be very helpful.

While this solution is nice and consistent our design problem is with the distinction of MPC vs LagrangeMultipliers vs Penalty. Conceptually the 3 of them are "constraints" and the normal user would understand them to belong to the ConstraintList. However in the practice the behaviour and construction of Master-Slave MPCs and of other constraints is completely different. Indeed the Elements/Condition interface is already ok for LagrangeMultipliers and PenaltyConstraints.

As you say the though Element/Condition is already ok for LM and Penality type of approaches, there is nothing restricting in applying them or formulating them through the ContraintEquations. Please correct me if I am wrong here.

If there is a possibility of applying LM and Penality with through ContraintEquations this would solve the confusion of how to treat LM & Penality & Master-Slave.

If there is a possibility, will this be a solution which can be considered ?

[RiccardoRossi]

Hi @adityaghantasala the comment you are making about the possibility of having LM and Penalty in the Constraints, is indeed also @pooyan-dadvand 's idea.

My point is that

1. MPCs are strictly linear, while LM and penalty are not. This means that LM/penalty is needed for any type of nonlinerar constraint (the simplest i can think of, is a is a large rotation)
2. LM and penalty lead to an "element" contribution that follows the standard FEM procedure. Conversely MPCs can be applied either by modifying the assembly phase (as you did) or "a posteriori" as a sparse-matrix-matrix product. Both possibilities need to remain open.

3. Ultimately there is one single type of MPC , providing a (linear) relation of the type

   y = T*x + b

   where x is input and y is output. Different physics are handled by providing different values for T and b (as well as for x and y). This is fundamentally different from elements...

4. Imagine now MPCs and LMs are in the same list. in the Build you would have to loop through that list, then when you encounter an LM, at the moment of assembling it you would have to take into account the MPCs it is subjected to (of course here you would have had to do previously a similar loop to mount auxiliary data structures for MPCs. I would say that all the looping would become very awkward.

One option indeed would be to have Elements, Conditions, MPCs, Constraints with the latter behaving exactly like an Element and condition. While this can be done ... the bookkeeping effort in the modelpart keeps growing, for something that is nothing else of an element in disguise...

[loumalouomega]

@RiccardoRossi just to say that I think you can compute penalty and LM without implement on elements, but when assembling the system and computing the constraints.

[RiccardoRossi]

@loumalouomega of course by functional programming you could achieve that. However let's just imagine those 3 scenarios:

• you want one point to move at a constant distance from another
• you want to apply a contact condition (jeje, teaser)
• you want to apply a minimal constraint for an RVE.

all 3 cases are perfectly possible with LM (but not with MPC). How would you instruct the builder that you want to do the first, but not the second? (or any combination of them for the matter)

[loumalouomega]

Ok, I got the idea

[RiccardoRossi]

@adityaghantasala do you agree with my post?

[RiccardoRossi]

I would make this change once the pybind is merged and the waves are settled

[edit:] i posted in the wrong windows...forget about this post

[adityaghantasala]

@RiccardoRossi I see the point. Just something that came to my mind : In the current setup , in the case of having two constraints, one LM and one penalty how would we specify the order?

Also I was thinking why not having linearized MPC, that is what if the constraint implementation provides the linearized matrices for the constraints ? then it will be possible to have a "non-linear" MPC right ?

[RiccardoRossi]

@adityaghantasala if you put LM as "element" than the ordering is guaranteed. If you mix it with MPCs you have an ordering problem

regarding your second point, think of a finite rotation. You do not want to write it in terms of incremental small rotations otherwise errors sum up. It in instead much better to have a residual and to write the problem as a residual minimization iteration, since at convergence the nonlinear error is going to zero, which is not guaranteed for the incremental form. This is actually the reason for which Total Lagrangian (or Update Lagrangian) formulations are superior to corotational ones.

[RiccardoRossi]

Hi guys, we have been discussing in the @KratosMultiphysics/technical-committee about adding MasterSlaveConstraints to the core.

We reached the agreement that 1 - a "MasterSlaveConstraint" should be a class, with an associated Id and flags (like elements and conditions but with different interface). The interface must allow somone to implement for example a non-linear link (to make an example a link active only in compression but not in traction). It will be however DIFFERENT from Elements and Conditions, implying that if someone will want to implement a LagrangeMultiplier constraint or a penalty constraint it will have to put it in the Elements list) 2 - It should be stored in the modelpart and submodelparts equivalently to Elements and Conditions. This means that once something is added to a submodelpart, it must also be automatically added to the parent modelparts, as happens for elements and conditions 3 - an interface equivalent to Elements and Conditions should be added to the modelpart

@adityaghantasala @philbucher can i count with your help on this? this would have to go in as a separate PR prior to the MPC branch being merged

[RiccardoRossi]

No,

I am speaking of the MPCs

On Sat, Apr 7, 2018 at 12:09 PM, Vicente Mataix Ferrándiz < notifications@github.com> wrote:

@RiccardoRossi https://github.com/RiccardoRossi I have a base class for all my conditions called PairedCondition

https://github.com/KratosMultiphysics/Kratos/blob/master/applications/ ContactStructuralMechanicsApplication/custom_conditions/paired_condition.h

It stores the pointer to the paired geometry, this is what do you mean?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/KratosMultiphysics/Kratos/issues/1698#issuecomment-379458602, or mute the thread https://github.com/notifications/unsubscribe-auth/AHr7EbYchN_cVbIISC7T412cGbFhEOZWks5tmJBSgaJpZM4Suv3c .

--

Riccardo Rossi

PhD, Civil Engineer

member of the Kratos Team: www.cimne.com/kratos

Tenure Track Lecturer at Universitat Politècnica de Catalunya, BarcelonaTech (UPC)

Full Research Professor at International Center for Numerical Methods in Engineering (CIMNE)

C/ Gran Capità, s/n, Campus Nord UPC, Ed. B0, Despatx 102

(please deliver post to the porters of building C1)

08034 – Barcelona – Spain – www.cimne.com -

T.(+34) 93 401 56 96 skype: rougered4

http://www.cimne.com/

https://www.facebook.com/cimne http://blog.cimne.com/ http://vimeo.com/cimne http://www.youtube.com/user/CIMNEvideos http://www.linkedin.com/company/cimne https://twitter.com/cimne

Les dades personals contingudes en aquest missatge són tractades amb la finalitat de mantenir el contacte professional entre CIMNE i voste. Podra exercir els drets d'accés, rectificació, cancel·lació i oposició, dirigint-se a cimne@cimne.upc.edu. La utilització de la seva adreça de correu electronic per part de CIMNE queda subjecte a les disposicions de la Llei 34/2002, de Serveis de la Societat de la Informació i el Comerç Electronic.

Imprimiu aquest missatge, només si és estrictament necessari. http://www.cimne.com/

[loumalouomega]

Yes, I thought about that, and I realize that my comment was meaningless, that why I removed, sorry for the noise

[adityaghantasala]

@RiccardoRossi Thanks for taking it up in the committee. Of course count me in.

While I could follow the part of adding the MasterSlaveConstraint to the ModelPart and the interface, few things which are not yet clear are :

• How do we store the MasterSlave Relation ? Also this was one of the main concern that was expressed in the comments above. Was there also a discussion about this ?

• Are the builder and solvers going to be merged into one ?

• Once the above two are clear, I can add a draft of this interface in a separate PR and we can further iterate/discuss.

• You mentioned that MPC is strictly linear in one of your comments above. In your latest comments you mentioned that one should be able to implement a non-linear link. I am a bit confused.

• Also I read a paper where there were using a non-linear MasterSlave relation and implement it using MPC by making use of the tangent matrix of the constraint relation. I still have to try it myself, but it seems it is possible to do.

[RiccardoRossi]

@adityaghantasala

The "Constraint" object will have to contain the relation between master and slave. It will be then used to construct the datastructures you need. Essentially what i am telling is that your constraint container should go in the modelpart and provide an hierarchical interface similar to elements. If you agree we can have a skype call to clarify this.

The long term goal is to have everything in a single builder and solver (well, a few of them, but not versions with and without MPCs). Having said this we will start with two separate versions and once we are sure everything works in production we'll preseve just one.

About the linearity and nonlinearity the @KratosMultiphysics/technical-committee pointed out during the meeting that i was wrong. For example one can implement a non-holonomic constraint (one directional constraint) by simply switching on and off the coefficients depending on wether compression or tension is found in the link. That is already an example of nonlinear constraint...which tells that i was wrong in my statement

[RiccardoRossi]

As a first step towards the implementation, the methods to be implemented for the python interface will be

    .add_property("MultiPointConstraints", ModelPartGetMultiPointConstraints1, ModelPartSetMultiPointConstraints1)
    .def("GetMultiPointConstraint", ModelPartGetMultiPointConstraint1)
    .def("GetMultiPointConstraint", ModelPartGetMultiPointConstraint2)
    .def("GetMultiPointConstraints", ModelPartGetMultiPointConstraints1)
    .def("SetMultiPointConstraints", ModelPartSetMultiPointConstraints1)
    .def("GetMultiPointConstraints", ModelPartGetMultiPointConstraints2)
    .def("SetMultiPointConstraints", ModelPartSetMultiPointConstraints2)
    .def("RemoveMultiPointConstraint", ModelPartRemoveMultiPointConstraint1)
    .def("RemoveMultiPointConstraint", ModelPartRemoveMultiPointConstraint2)
    .def("RemoveMultiPointConstraint", ModelPartRemoveMultiPointConstraint3)
    .def("RemoveMultiPointConstraint", ModelPartRemoveMultiPointConstraint4)
            .def("RemoveMultiPointConstraints", &ModelPart::RemoveMultiPointConstraints)
    .def("RemoveMultiPointConstraintFromAllLevels", ModelPartRemoveMultiPointConstraintFromAllLevels1)
    .def("RemoveMultiPointConstraintFromAllLevels", ModelPartRemoveMultiPointConstraintFromAllLevels2)
    .def("RemoveMultiPointConstraintFromAllLevels", ModelPartRemoveMultiPointConstraintFromAllLevels3)
    .def("RemoveMultiPointConstraintFromAllLevels", ModelPartRemoveMultiPointConstraintFromAllLevels4)
            .def("RemoveMultiPointConstraintsFromAllLevels", ModelPartRemoveMultiPointConstraintsFromAllLevels)
    .def("MultiPointConstraintsArray", &ModelPart::MultiPointConstraintsArray, return_internal_reference<>())
    .def("CreateNewMultiPointConstraint", ModelPartCreateNewMultiPointConstraint)
            .def("AddMultiPointConstraint", &ModelPart::AddMultiPointConstraint)
            .def("AddMultiPointConstraints",AddMultiPointConstraintsByIds)

these are copied from the Element interface of the modelpart.

A crucial aspect is that all MPCs will have to be created by the CreateNewMultiPointConstraint function, of which we'll have to define multiple versions.

What i am thinking of is to take as inputs

• slave_dofs
• master_dofs
• CoefficientMatrix (of size slave_dofs.size() * master_dofs.size()
• CoefficientVector (of size slave_dofs.size() )

But i guess there will be variations.

Can u explain what is the current interface?

[adityaghantasala]

@RiccardoRossi Can we discuss this further tomorrow during the regular skype meetings ?? Then I can start working on the interface ?

[RiccardoRossi]

@adityaghantasala unfortunately tomorrow i have another meeting at 11. if you like we can make it at 10.

[adityaghantasala]

Yes we can do it at 10. @philbucher is this also ok for you ?

[philbucher]

Fine for me and also for the others I think. @loumalouomega FYI

[loumalouomega]

Yep

El lun., 9 abr. 2018 a las 17:09, Philipp Bucher (notifications@github.com) escribió:

Fine for me and also for the others I think. @loumalouomega https://github.com/loumalouomega FYI

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/KratosMultiphysics/Kratos/issues/1698#issuecomment-379785790, or mute the thread https://github.com/notifications/unsubscribe-auth/ATEMgA17bVi6tFaHjLJ4x-TbwNJlNGa6ks5tm3lxgaJpZM4Suv3c .

[philbucher]

We are ready

[RiccardoRossi]

ok, this has been discussed and decided. closing

[adityaghantasala]

I am reopening this to discuss more about the interface of the constraints in the modelpart. This way all the discussion is at one place. I was looking at necessary functions at the interface. I have these following points/questions

• There should be some additional functions (in addition to the ones posted above).
• These are mainly the functions to add / remove / modify the masters and their weights
• Setting the constraint active / inactive is also important. This can be done via the Transformation matrix it returns but it will save a few operations if there is a ACTIVE flag on the constraint (like on an element/condition)
• There will be ConstraintsContainerType like the ElementsContainerType (correct me if I am wrong). Many times to build a proper matrix for a constraint, it is necessary if a equation Id related to a DOF is a slave or not. So is it possible to overwrite the = operator for the constraint base and use it to compare both the slave dof and its equationID (in builder and solver its easy to work with equation ids )
• Removing the constraints should always be done from all the levels. This is necessary to avoid any leftover constraints on some levels effecting the solution.
• applying the constraint T matrix (s = T*m+c) to the element matrices or condition matrices should be done in a small function inside the builder and solver. Is this ok ?

[RiccardoRossi]

Hi aditya, i'll try to answer point by point

1. I guess we will need some additional functions, which match the Elements and Conditions interface, however i guess that the functions you want will be functions of the constraint itself rather than of the modelpart. We can derive the MPC object from flags so that it will have the flags interface
2. I don't follow what you are asking about EquationIds. Within the Build() function it is ok to use EquationIds, however to mount the graph etc you will also need the dofs. (Node and variable is a way of specifying one). Can u specify a bit more what is needed?
3. about removing, it would be nice to have the capability of calling a strategy with a modelpart that is not the root, for example to solve just one portion of the domain. Why do you foresee problems with this?
4. yes, the builder and solver will ask T and c to the list of active constraints and use it for mounting the needed data structures (for example IT COULD assemble a sparse matrix Tsparse as an alternative to the elemental formulation. Note that i am NOT telling that you have to do it that way, simply that we have alternatives in how the builder and solver operates. ) The final application will be then done in the builder and solver

[adityaghantasala]

Sorry for late reply.

I guess we will need some additional functions, which match the Elements and Conditions interface, however i guess that the functions you want will be functions of the constraint itself rather than of the modelpart. We can derive the MPC object from flags so that it will have the flags interface

agreed. I was also thinking like this.

I don't follow what you are asking about EquationIds. Within the Build() function it is ok to use EquationIds, however to mount the graph etc you will also need the dofs. (Node and variable is a way of specifying one). Can u specify a bit more what is needed?

Looking into the scenario of element wise operation, once we get the element matrices from any particular element, we should know if we need to change (we need to change only if a node of this element is marked slave) these matrices to apply any MPC to it. Here there is a check (if node is slave) and accessing the container of MasterSlaveConstraints to get the T and c of that particular constraint of that node. Here since equation Ids are already available it would be easy to access the container of MasterSlaveConstraints via equation ID.

about removing, it would be nice to have the capability of calling a strategy with a modelpart that is not the root, for example to solve just one portion of the domain. Why do you foresee problems with this?

A bit diffused but yes. Lets say we remove the constraint only from a submodelpart this means the modelpart which is used by strategy will still have it and builder and solver does not see that this constraint is removed and will apply it anyway.

yes, the builder and solver will ask T and c to the list of active constraints and use it for mounting the needed data structures (for example IT COULD assemble a sparse matrix Tsparse as an alternative to the elemental formulation. Note that i am NOT telling that you have to do it that way, simply that we have alternatives in how the builder and solver operates. ) The final application will be then done in the builder and solver

agreed, this is also what I have in mind.

[RiccardoRossi]

Hi @adityaghantasala

essentially i foresee having 2 stages in the builder and solver: stage1 - BuilderAndSolver queries the MPCs for their coefficients Matrix and Vector and uses it for building an auxiliary data structure (owned solely by the B&S) stage2 - During assembly the builder and solver will query its auxiliary data structure and do all the necessary ops (or do nothing as needed)

regarding the removal or addition of MPCs this is similar to the addition or removal of elements. I mean, if you change something that affects the graph, you'll have to call the method Clear of the strategy for it to work correctly....this is a user-side thing as it is for elements

[adityaghantasala]

regarding the removal or addition of MPCs this is similar to the addition or removal of elements. I mean, if you change something that affects the graph, you'll have to call the method Clear of the strategy for it to work correctly....this is a user-side thing as it is for elements

If we want to do it like that, I have no problem. I was thinking of being more careful.

essentially i foresee having 2 stages in the builder and solver: stage1 - BuilderAndSolver queries the MPCs for their coefficients Matrix and Vector and uses it for building an auxiliary data structure (owned solely by the B&S) stage2 - During assembly the builder and solver will query its auxiliary data structure and do all the necessary ops (or do nothing as needed)

I see the two stages. In fact, in my previous implementation there are these two stages, and I was storing this data structure in the constraint itself. I am now thinking if it is effective to combine these two stages and have it in the one stage during the assembly.

In the mean time I am also preparing a draft implementation of the interface based on the discussions we had and inputs from here.

[RiccardoRossi]

@adityaghantasala i think that the bottomline is that MPCs in modelpart have the task of storing the information and making it available to the B&S. The latter is in charge of making use of such information.

I guess the easiest is if we have a video call to speak of the details

[adityaghantasala]

Hi @RiccardoRossi and @pooyan-dadvand

I implemented the suggested and discussed interface of master-slave constraints to ModelPart. I made a pull request to make the changes easily visible.

I retained the previous implementation of storing the constraint (using key and ID) for time being.

The pull request is here #2204

[philbucher]

I am closing since the coresponding pr is merged

Please repopen if I missed sth