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 <[email protected]> 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: >>>>> [email protected] >>>>> 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/fo >>>>> rum/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 [email protected]. >>>>> 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/fo >>> rum/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 [email protected]. >>> 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 [email protected]. > 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 [email protected]. For more options, visit https://groups.google.com/d/optout.
