Cantilever Beam...more complicated than expected

[Nicholas Padon]

I'm playing with dealii (8.3.0) and trying to solve a 2-D, plane strain elasticity problem for a simple cantilever beam (say 5x1). Although I acknowledge that dealii is made for more sophisticated problems than a "stupid" beam, I think that the use (and documentation) of such simple problems might increase dealii's broader appeal (for example to dumb engineers like myself).

Although I follow step-8 and get expected results for distributed forces (i.e. right hand sides that are continuous functions), I'm having a little difficulty finding a clean way to represent nodal forces (i.e. point forces). My RightHandSide function looks like:

void RightHandSide<dim>::vector_value(const Point<dim>& p,

Vector<double>& values) const

{

Assert(values.size()==dim,ExcDimensionMismatch(values.size(),dim));

Assert(dim ==2,ExcNotImplemented());

//apply a force at the top end of the beam

if(p(1)==1 && p(0)==5)

{

values(0)=0;

values(1)=-1000;

}

else

{

values(0)=0;

values(1)=0;

}

}

When iterating through each cell using Gaussian quadrature to assemble the stiffness matrix and right hand side (following step-8), any point force (say, one applied directly at a node for simplicity) isn't "seen" by the quadrature points, so the right hand side ends up being zero. Fine. So I've got to update the nodal forces some other way...

Let's say I've got a list of vertices (say from triangulation.get_vertices()) and I can return the associated force at each node via RightHandSide.vector_value(). My problem is figuring out how to associate the vertex/point information (and the forces associated with that vertex) with the appropriate global degree of freedom such that I can get an updated right hand side.

I've (partially) tried a couple of things:

1. Using DoFTools::map_support_points_to_dofs in some way, but, as discussed in a previous google group post, this function should probably be deprecated (or updated to return a multimap) and thus is currently not a good candidate for choosing the nodal forces.

2. Iterating through active cells, adding each cell's vertex information (Point<dim>) and global dof information to a unordered_map of type <Point<dim>,std::vector<types::global_dof_index>>. Then, at the end of assemble_system(), I write the value of RightHandSide at each point in the unordered_map to the location specified by the point's global dofs in the system_rhs vector. This is kind of a dumb method (and it requires me to implement a hash function and point comparison function to make the unordered_map work) but it seems to function... 

3. Using VectorTools::create_point_source_vector (instead of my RightHandSide function) but I get an error that the function is only implemented in 1-D.

4. Using VectorTools::create_right_hand_side with a custom quadrature rule and mapping (i.e. quadrature points are actually just the vertices, with weights equal to one for each points - and the mapping returns just the function values). Note that I haven't implemented this option because it just seems too crazy to do such a simple thing...

My question is: is there a more elegant way to represent point forces than what I've described in Option 2?

DealII is pretty great, by the way.

P.S. I attach my current source code, although I acknowledge that I haven't yet updated the cell_matrix to reflect the plane strain stiffness matrix (right now it's the perfect 2-D solid, as Wolfgang noted in a previous tutorial). Thanks in advance.

[Bruno Turcksin]

Hi,

On Sunday, October 4, 2015 at 9:08:32 PM UTC-5, Nicholas Padon wrote:

I'm playing with dealii (8.3.0) and trying to solve a 2-D, plane strain elasticity problem for a simple cantilever beam (say 5x1). Although I acknowledge that dealii is made for more sophisticated problems than a "stupid" beam, I think that the use (and documentation) of such simple problems might increase dealii's broader appeal (for example to dumb engineers like myself).

Although I follow step-8 and get expected results for distributed forces (i.e. right hand sides that are continuous functions), I'm having a little difficulty finding a clean way to represent nodal forces (i.e. point forces). My RightHandSide function looks like:

void RightHandSide<dim>::vector_value(const Point<dim>& p,

Vector<double>& values) const

{

Assert(values.size()==dim,ExcDimensionMismatch(values.size(),dim));

Assert(dim ==2,ExcNotImplemented());

//apply a force at the top end of the beam

if(p(1)==1 && p(0)==5)

 

This cannot work. There is almost no chance that p(1) will be exactly 1, it will probably something like 1.00000000001 so you need to do something like this if((std::abs(p(1)-1.)<esp) && std::abs(p(1)-5.)<eps) with eps small (for ex. 10^{-12})

My question is: is there a more elegant way to represent point forces than what I've described in Option 2?

If I understand correctly, your right-hand-side is similar to a Dirac delta. So when you write your weak form, your right-hand-side is something like \int \int \phi(x,y) \delta(x_0,y_0) dx dy = \phi(x_0,y_0) So to build your right-hand-side, you just need to evaluate the test functions at one point.

Best,

Bruno

[Nicholas Padon]

If we ignore the floating point comparison problem, and focus on your suggestion: when iterating through the cells, I would thus have to include the cell vertices as evaluation points (similar to my option 4), and I would have to flag the vertices as "touched" in some way such that I don't count a particular vertex twice (and thus include a nodal force twice). Is that "better" than Option 2 (i.e. directly updating the system_rhs)? 

[Bruno Turcksin]

What ever works is the best option. Now, if you use Lagrange finite elements (FE_Q), then you know that at each vertex there is only one non-zero basis function and the value is one. Now using the numbering of FE_Q (https://dealii.org/8.3.0/doxygen/deal.II/classFE__Q.html), you can do something like:

for (loop on every cell)
    for (loop on every vertex of the cell)
        if (vertex is in the list of discrete points) // here you need to implement a point comparison be careful about floating point comparison
        {
            cell->get_dof_indices(global_dof_indices);
            // Let's assume that you are on 2D and the top left vertex is the one you are interested in
            rhs[global_dof_indices[2]] = value_at_this_vertex; // You set right the value, not add it so there is no need to know this vertex has already be done
        }

This is basically the same as your option 2 but there is not need for an unordered map.

Best,

Bruno

--
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 a topic in the Google Groups "deal.II User Group" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dealii/HTcSxhnzc5s/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dealii+un...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Nicholas Padon]

That does look like a better option (particularly as it does away with the need for an unordered map). Thanks for your help.

Nick

[Nicholas Padon]

I include what I finally did, which is:

...loop on cells

// Assembling the cell rhs - nodal values first

            for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell;++v)

            {

            double global_vertex_location = cell->vertex_index(v);

            AssertIndexRange(global_vertex_location,vertices_touched.size());

            if(vertices_touched[global_vertex_location]==false)

            {

vertices_touched[global_vertex_location] = true;

right_hand_side.vector_value(cell->vertex(v),rhs_nodal_values);

for(unsigned int i=0;i<fe.dofs_per_vertex;++i)

cell_rhs(fe.dofs_per_vertex*v+i) = rhs_nodal_values(i);

            }

            }

It actually corresponds to the procedure indicated in the glossary entry for user flags, where I have a vector<bool> for each vertex used in the triangulation. In this way, I don't have to have a point comparison function at all, which seems even better. Thanks.

[Jiaming Kang]

Very good example. I am learning dealii to do some simple structural analysis as well. It is quit useful for me. Thanks.