See step-27 <https://www.dealii.org/developer/doxygen/deal.II/step_27.html>:
VectorTools::interpolate_boundary_values (dof_handler,
> 0,
> ZeroFunction<dim>(),
> constraints);
On Friday, October 7, 2016 at 2:02:18 PM UTC+2, Deepak Gupta wrote:
>
> Thanks JP,
>
> I took the above points into my code except that I do not know/understand
> how to add boundary conditions in the constraint matrix without using the
> apply_boundary_values. Could you divert me to an example or discussion
> where this is done.
>
> Best
> Deepak
>
> On Fri, Oct 7, 2016 at 1:33 PM, Jean-Paul Pelteret <jppelte...@gmail.com>
> wrote:
>
>> Hi Deepak,
>>
>> Its not pretty but at a glance this does look plausibly correct. There
>> are a few more things that you should check:
>> - Add the boundary conditions directly to the constraints matrix, instead
>> of using MatrixTools::apply_boundary_values. You can just use a
>> ZeroFunction on boundary 42 instead of the function you've defined here.
>> - You've created a temporary quadrature object "QGauss<2>
>> temp_quad(cell->active_fe_index()+3);" which is suspicious. Remove this.
>> You can retrieve the active quadrature formula from FEValues.
>> - Run this in debug mode to check that there are no other mistakes that
>> we can't pick up by just looking at the code.
>>
>> If you do all of this and still have no success, then post the updated
>> code and I'll fiddle with it later.
>> J-P
>>
>>
>> On Friday, October 7, 2016 at 12:19:39 PM UTC+2, Deepak Gupta wrote:
>>>
>>> Thanks JP for the elaborate clarification. I quickly tried to test a
>>> *distributed
>>> source term* commenting out the point load function and adding the
>>> following before distribute_local_to_global:
>>>
>>> //Calculating cell_rhs
>>> for(unsigned int i = 0; i < dofs_per_cell; ++i){
>>> const unsigned int component_i =
>>> cell->get_fe().system_to_component_index(i).first;
>>> for(unsigned int q_point = 0; q_point < n_q_points;
>>> ++q_point){
>>> double body_load = 0;
>>> if (component_i == 1 && cell_itr%4 == 3) body_load = 0.5;
>>> cell_rhs(i) += fe_values.shape_value(i, q_point) *
>>> body_load * //rhs_values[q_point](component_i) *
>>> fe_values.JxW(q_point);
>>> }
>>> }
>>>
>>>
>>> However, the results are still similar, although magnitudes differ
>>> (globally_p_2.png and locally_p_3.png). The cases with global and local
>>> refinement are attached. The strange thing is that the symmetry in absolute
>>> value of x-displacement is lost when one element is of higher order. I
>>> would expect a bit of difference but not this much.
>>>
>>> Do you think the source term implementation is correct? I added force
>>> term only in the two cells on the ride side.
>>>
>>> Best regards
>>> Deepak
>>>
>>> On Fri, Oct 7, 2016 at 11:48 AM, Jean-Paul Pelteret <
>>> jppelte...@gmail.com> wrote:
>>>
>>>> Ok, so its a benchmark test not one you've come up with yourself. Fair
>>>> enough...
>>>>
>>>> What I'm suggesting is that you implement the RHS contribution from the
>>>> load within the local assembly loop. You can have a look at how the point
>>>> source is implemented in VectorTools::create_point_source_vector and
>>>> duplicate this within your assembly routine (you could probably either do
>>>> this localised to one cell or distribute the load contribution over the
>>>> neighbouring cells). This way when you call distribute_local_to_global
>>>> this
>>>> contribution gets scaled correctly. Why you can't do this with your
>>>> current
>>>> function "add_point_source_to_rhs" is because it directly modifies the
>>>> global system RHS vector.
>>>>
>>>> I hope that this explains better what I think you might need to do to
>>>> correct the issue.
>>>>
>>>> J-P
>>>>
>>>> On Friday, October 7, 2016 at 11:26:18 AM UTC+2, Deepak Gupta wrote:
>>>>>
>>>>> Thanks JP. Indeed a concise example should help and that is where I
>>>>> started. I took a scalar-valued example, however, it did not show the
>>>>> same
>>>>> problem as this one. Thus, I decided to move towards this case. I take
>>>>> this
>>>>> blame that I should have an even simpler case, and I hope I can come up
>>>>> soon.
>>>>>
>>>>> In the meantime, I will take the suggestions into account and see if
>>>>> the issue goes away.
>>>>>
>>>>>
>>>>>
>>>>> *"Anyway, I've noticed that you call "add_point_source_to_rhs" after
>>>>> you've applied the hanging node constraints and boundary conditions via
>>>>> hanging_node_constraints.distribute_local_to_global. This might be the
>>>>> problem because the matrix and RHS contribution is modified when the
>>>>> local
>>>>> constraints are applied. I would suggest that you implement this as a
>>>>> standard traction boundary condition and then add its local contribution
>>>>> to
>>>>> the global RHS vector as per usual. This would eliminate another possible
>>>>> source of error."*
>>>>> So, would it be correct, if I move my add_point_source_to_rhs function
>>>>> above create_local_to_global? Also, to compare with the benchmarks, I
>>>>> need
>>>>> to keep it a point load. Would that still be possible with traction
>>>>> boundary condition?
>>>>>
>>>>> Best,
>>>>> Deepak
>>>>>
>>>>> On Fri, Oct 7, 2016 at 11:06 AM, Jean-Paul Pelteret <
>>>>> jppelte...@gmail.com> wrote:
>>>>>
>>>>>> Dear Deepak,
>>>>>>
>>>>>> Its always easier to help debug a well document, clear and concise
>>>>>> example. It was less about checking the implementation than
>>>>>> understanding
>>>>>> what the code is doing and how all of the parts interact. After all, we
>>>>>> can
>>>>>> only work with what you give us.
>>>>>>
>>>>>> Anyway, I've noticed that you call "add_point_source_to_rhs" after
>>>>>> you've applied the hanging node constraints and boundary conditions via
>>>>>> hanging_node_constraints.distribute_local_to_global. This might be
>>>>>> the problem because the matrix and RHS contribution is modified when the
>>>>>> local constraints are applied. I would suggest that you implement this
>>>>>> as a
>>>>>> standard traction boundary condition and then add its local contribution
>>>>>> to
>>>>>> the global RHS vector as per usual. This would eliminate another
>>>>>> possible
>>>>>> source of error.
>>>>>>
>>>>>> Yes, as best I can discern VectorTools::create_point_source_vector
>>>>>> should do the same job for you (i.e. create the equivalent RHS
>>>>>> contribution
>>>>>> as a point load). I've never used it though so I don't know what the
>>>>>> implications of using it in this context are (i.e. the interaction of
>>>>>> hp-constraints).
>>>>>>
>>>>>> Regards,
>>>>>> J-P
>>>>>>
>>>>>> On Friday, October 7, 2016 at 8:41:42 AM UTC+2, Deepak Gupta wrote:
>>>>>>>
>>>>>>> Dear JP,
>>>>>>>
>>>>>>> Thanks for looking at the example. The code itself is not a simple
>>>>>>> example that I implemented. Rather it is a part of my exact problem and
>>>>>>> I
>>>>>>> have cut it out from the rest of the code to demonstrate my problem. I
>>>>>>> am
>>>>>>> working on design optimization and this part goes in there as the
>>>>>>> modeling
>>>>>>> part.
>>>>>>>
>>>>>>> Indeed the compilation of the B-matrix might look complicated.
>>>>>>> However, a closer look might reveal that B matrix is compiled as shown
>>>>>>> in *http://homepages.rpi.edu/~des/4NodeQuad.pdf
>>>>>>> <http://homepages.rpi.edu/~des/4NodeQuad.pdf>*, page 2, slide 3.
>>>>>>> This link is just a random example I took. I hope there is no error in
>>>>>>>
>>>>>>> that part of the code since I have been using it and it worked in the
>>>>>>> past.
>>>>>>> Also, if that formulation would be wrong, I would wonder why the
>>>>>>> globally
>>>>>>> fine mesh would give the correct results. The values 420, 428 etc are
>>>>>>> correct since they are very close to a MATLAB version of FEM code that
>>>>>>> I
>>>>>>> compared with.
>>>>>>>
>>>>>>> Thanks for pointing out the *vertex_dof_index* function. However, I
>>>>>>> have some question related to it now:
>>>>>>>
>>>>>>> 1) When I describe the point load using my approach, I get the load
>>>>>>> where I wanted (see my_approach_point_load.png).
>>>>>>>
>>>>>>> 2) However, when I use *vertex_dof_index *function, the load gets
>>>>>>> attached to a different component of a wrong vertex (see
>>>>>>> new_point_load.png). I am wondering how I should use this function, or
>>>>>>> may
>>>>>>> be there is something else wrong with my code that makes it wrong.
>>>>>>> I used following lines in the code:
>>>>>>> system_rhs(cell->vertex_dof_index(v, 0)) = 0;
>>>>>>> system_rhs(cell->vertex_dof_index(v, 1)) = 1;
>>>>>>>
>>>>>>> 3) In addition, I would like to know whether I can use
>>>>>>> *Ve*
>>>>>>> *ctorTools::create_point_source_vector *for the same? This function
>>>>>>> used to help me for bilinear elements.
>>>>>>>
>>>>>>> Best regards
>>>>>>> Deepak
>>>>>>>
>>>>>>> On Thu, Oct 6, 2016 at 8:10 PM, Jean-Paul Pelteret <> wrote:
>>>>>>>
>>>>>>>> Hi Deepak,
>>>>>>>>
>>>>>>>> You've chosen a difficult example for us to evaluate, because your
>>>>>>>> implementation of the plane stress condition is rather hard to follow.
>>>>>>>> Can
>>>>>>>> you provide a reference from which you're implementing this?
>>>>>>>>
>>>>>>>> At the very least your implementation of "add_point_source_to_rhs"
>>>>>>>> looks suspicious to me. The ordering of the vertices has, in general,
>>>>>>>> nothing to do with the ordering of the DoFs. You should use the
>>>>>>>> vertex_dof_index
>>>>>>>> <https://www.dealii.org/developer/doxygen/deal.II/classDoFAccessor.html#a5560151b5407e4851d5c1009c7753764>
>>>>>>>> function
>>>>>>>> to retrieve the DoF associated with a vertex. There are a number
>>>>>>>> of posts
>>>>>>>> <https://groups.google.com/forum/#!searchin/dealii/vertex_dof_index%7Csort:relevance>
>>>>>>>>
>>>>>>>> that discuss the use of this function.
>>>>>>>>
>>>>>>>> J-P
>>>>>>>>
>>>>>>>> On Thursday, October 6, 2016 at 7:39:28 PM UTC+2, Deepak Gupta
>>>>>>>> wrote:
>>>>>>>>>
>>>>>>>>> Hi Wolfgang,
>>>>>>>>>
>>>>>>>>> Yes it is ConstraintMatrix::distribute_local_to_global.
>>>>>>>>>
>>>>>>>>> I have finally come up with a small example that clearly shows the
>>>>>>>>> issue I am facing. Attached is the code and 3 outputs. I have a
>>>>>>>>> rectangle
>>>>>>>>> domain of 4X8 elements. A point load (0, 1) is applied at the second
>>>>>>>>> vertex
>>>>>>>>> from the top at the right boundary and the left boundary is fixed.
>>>>>>>>> The
>>>>>>>>> y-displacement will be maximum at the load point. Now I look at 3
>>>>>>>>> cases:
>>>>>>>>>
>>>>>>>>> 1) Using *globally biquadratic elements* -> displacement value is
>>>>>>>>> *420.2
>>>>>>>>> *(as in biquadratic.png)
>>>>>>>>> *.*2) Using* globally bicubic elements* -> displacement value is
>>>>>>>>> *426.8* (as in bicubic.png).
>>>>>>>>> 3) Using *biquadratic for all elements*, except that the *element
>>>>>>>>> with index 6 is bicubic*, displacement value is *61.45 *(as in
>>>>>>>>> mix.png).
>>>>>>>>>
>>>>>>>>> I would expect the result of *3) *to be between *1) *and *2)
>>>>>>>>> *ideally.
>>>>>>>>> I have tested it for finer meshes, but such issues exists. Could
>>>>>>>>> someone
>>>>>>>>> help me in locating the problem. Thanks in advance.
>>>>>>>>>
>>>>>>>>> Best regards
>>>>>>>>> Deepak
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> On Thu, Oct 6, 2016 at 7:16 PM, Wolfgang Bangerth <> wrote:
>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> Deepak,
>>>>>>>>>>
>>>>>>>>>> By the way, I am constructing the cell matrices using B'DB, where
>>>>>>>>>>> B is
>>>>>>>>>>> the strain-displacement matrix and D is the constitutive
>>>>>>>>>>> relation. Then
>>>>>>>>>>> I add it to the system matrix using distribute_local_to_gobal
>>>>>>>>>>> function.
>>>>>>>>>>> Since B and D are the same for global and adaptive case, I do
>>>>>>>>>>> not expect
>>>>>>>>>>> any difference in the cell matrices. The system_rhs is correct
>>>>>>>>>>> since I
>>>>>>>>>>> have only a point load and I can check it by printing the
>>>>>>>>>>> vector. Thus,
>>>>>>>>>>> the possible issue can be only in the assembly (for the local
>>>>>>>>>>> refinement
>>>>>>>>>>> one). My problem is vector-values. In case this information can
>>>>>>>>>>> lead to
>>>>>>>>>>> any other suggestion from someone, it would be great.
>>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> I don't have any experience with this approach, but at the end of
>>>>>>>>>> the day, you get a local matrix A_K = B_K' D_K B_K for each cell
>>>>>>>>>> K, and
>>>>>>>>>> that should be all that matters. I assume you mean
>>>>>>>>>> ConstraintMatrix::distribute_local_to_gobal, not
>>>>>>>>>> cell->distribute_local_to_gobal, right?
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> Meanwhile I am going to extend the simple example for a basic
>>>>>>>>>>> version of
>>>>>>>>>>> the elastic problem I am trying to solve. Hope then I can figure
>>>>>>>>>>> out the
>>>>>>>>>>> error.
>>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> Good plan.
>>>>>>>>>>
>>>>>>>>>> Best
>>>>>>>>>> W.
>>>>>>>>>>
>>>>>>>>>> --
>>>>>>>>>>
>>>>>>>>>> ------------------------------------------------------------------------
>>>>>>>>>> Wolfgang Bangerth email:
>>>>>>>>>> bange...@colostate.edu
>>>>>>>>>> www:
>>>>>>>>>> http://www.math.colostate.edu/~bangerth/
>>>>>>>>>>
>>>>>>>>>> --
>>>>>>>>>> The deal.II project is located at http://www.dealii.org/
>>>>>>>>>> For mailing list/forum options, see
>>>>>>>>>> https://groups.google.com/d/forum/dealii?hl=en
>>>>>>>>>> --- You received this message because you are subscribed to the
>>>>>>>>>> Google Groups "deal.II User Group" group.
>>>>>>>>>> To unsubscribe from this group and stop receiving emails from it,
>>>>>>>>>> send an email to dealii+unsubscr...@googlegroups.com.
>>>>>>>>>> For more options, visit https://groups.google.com/d/optout.
>>>>>>>>>>
>>>>>>>>>
>>>>>>>>> --
>>>>>>>> The deal.II project is located at http://www.dealii.org/
>>>>>>>> For mailing list/forum options, see
>>>>>>>> https://groups.google.com/d/forum/dealii?hl=en
>>>>>>>> ---
>>>>>>>> You received this message because you are subscribed to the Google
>>>>>>>> Groups "deal.II User Group" group.
>>>>>>>> To unsubscribe from this group and stop receiving emails from it,
>>>>>>>> send an email to dealii+unsubscr...@googlegroups.com.
>>>>>>>> For more options, visit https://groups.google.com/d/optout.
>>>>>>>>
>>>>>>>
>>>>>>> --
>>>>>> The deal.II project is located at http://www.dealii.org/
>>>>>> For mailing list/forum options, see
>>>>>> https://groups.google.com/d/forum/dealii?hl=en
>>>>>> ---
>>>>>> You received this message because you are subscribed to the Google
>>>>>> Groups "deal.II User Group" group.
>>>>>> To unsubscribe from this group and stop receiving emails from it,
>>>>>> send an email to dealii+unsubscr...@googlegroups.com.
>>>>>> For more options, visit https://groups.google.com/d/optout.
>>>>>>
>>>>>
>>>>> --
>>>> The deal.II project is located at http://www.dealii.org/
>>>> For mailing list/forum options, see
>>>> https://groups.google.com/d/forum/dealii?hl=en
>>>> ---
>>>> You received this message because you are subscribed to the Google
>>>> Groups "deal.II User Group" group.
>>>> To unsubscribe from this group and stop receiving emails from it, send
>>>> an email to dealii+unsubscr...@googlegroups.com.
>>>> For more options, visit https://groups.google.com/d/optout.
>>>>
>>>
>>> --
>> The deal.II project is located at http://www.dealii.org/
>> For mailing list/forum options, see
>> https://groups.google.com/d/forum/dealii?hl=en
>> ---
>> You received this message because you are subscribed to the Google Groups
>> "deal.II User Group" group.
>> To unsubscribe from this group and stop receiving emails from it, send an
>> email to dealii+unsubscr...@googlegroups.com.
>> For more options, visit https://groups.google.com/d/optout.
>>
>
>
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
