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/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. >> > > -- 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.
