Re: [deal.II] Face and line orientation dependent dof permutations in FE_Q
On 1/13/21 12:32 PM, Konrad Simon wrote: I have a quick question: On complicated meshes some faces can have a non-standard orientation or they can be rotated against the neighboring face. This causes inconsistencies for some vector elements. In principle this potentially causes inconsistencies for FE_Q elements as well (for all dofs located at vertices, lines and faces). But that is being taken care of in the library, so everything seems fine. In order to understand the issue better for other elements: where in the library is this being taken care of? I have difficulties to find that. Remark: I found a table in (fe.h/fe.cc) that derived elements need to implement that takes care of local dof permutations on quads (not full faces and hence not line dofs or vertex dofs). Now, lines can have different orientations but lines also permute within a cell upon face orientation/rotation changes. Can anyone point me to the place in Deal.II where this is being taken care of for FE_Q? Konrad -- my apologies not only for not getting to this question quicker, but in fact for not helping you along over the past few months. I saw that you are working on this, but it's been a very long time since I understood these issues sufficiently well to help out and I haven't had the time to really read up on it again in enough detail to be useful. I do know that it's a topic for which we, collectively, have probably lost our expertise and that there is likely nobody else who knows it well either. You might be the expert on this at the moment. Fundamentally, the place to look is to search through the finite element classes where they ask for face and line orientations. It is possible that you don't actually have to do much there: Finite element classes generally work in the local (cell) coordinate system, where you enumerate shape functions as 0...dofs_per_cell. So when, for example, you compute the local matrix, it really is local. The interesting part happens when you translate what *global* degrees of freedom correspond to the local ones, i.e., when you call cell->get_dof_indices() (==DoFCellAccessor::get_dof_indices). This is where on every cell has to collect which DoFs live on it, and that has to take into account the orientations of lines and faces come into play. The cell DoF indices are cached, but you can see how this works in dof_accessor.templates.h, starting in line 902 using the process_dof_indices() functions. What this means is that if you have two cells that disagree about the relative orientation of their shared face, they when you call cell->get_dof_indices() on these two cells, you'd end up with a different ordering of the DoFs on the face. My usual approach to check this would be to create a mesh with two cells so that one has a non-standard face orientation, assign a Q3 element to the DoFHandler, and check that that is really what happens. (I generally turn these little programs into tests for the test suite, just so that what I just verified to be correct -- or not! -- is actually checked automatically in the future. I would very much love to see you submit such tests as well!) So this is the process for FE_Q, where there is really not very much to do, or in fact for all interpolatory elements. The situation becomes a bit more difficult for things such as the Nedelec element (also the Raviart-Thomas one) where the shape functions are tangential to the edge. Let's say you have the lowest order Nedelec element. Then there is just one DoF on each edge, and both adjacent cells agree about this. But they also have to agree *which direction the vector-valued shape function should point to*, and that's trickier. I can't remember how we do this (or whether we actually do this at all). I *think* that the RT element gets this right in 2d, but I don't remember about the 3d case. Again, it would probably be quite helpful to write little programs and turn them into tests. I know or strongly suspect that there are bugs for RT and Nedelec elements on meshes where faces have non-standard orientation. I've been thinking for years how I would go about debugging this, but I've never had the time to actually do it. But I'll share what I think would be the tools to do this: 1/ Creating meshes with non-standard orientations: You can feed Triangulation::create_triangulation() directly so that the faces have non-standard orientations with regard to each other. (See step-14 for how to call this function.) In 2d, this would work as follows: Think of a mesh like this: 3---4---5 | | | 0---1---2 You'd describe the two cells via their vertex indices as 0 1 3 4 1 2 4 5 but if you want a mismatched edge orientation, you'd just use 0 1 3 4 (thinks the edge is 1->4) 5 4 2 1 (thinks the edge is 4->1) or any other rotation. A good test would probably just try all 2x4 possible rotations of both cells. 2/ How do you check whether things are right
Re: [deal.II] Re: Transfer vector of solutions
HI Karthik,
Glad we could help :-)
To Question 1:
So you estimate the error for the second component of each of your
solutions, and then you add up the estimates on each cell. What is the
reasoning behind summing up the error of multiple solutions?
Assume that the error estimate of one of your solutions is larger than all
others by magnitudes, then this would be the dominant part of your estimate
sum, and the other solutions would have no influence at all.
Maybe it would be reasonable to only pick one of your estimates as a
measure for grid refinement.
First, you should find a detailed answer to the question why you want to
refine your grid, and then find a way how you want to do it.
To Question 2:
`GridRefinement::refine_and_coarsen_fixed_number` always refines and
coarsens a the specified fraction of ALL cells. With this you can control
the growth of the mesh size, rather than the reduction of the error. For
the latter, you can use `GridRefinement::refine_and_coarsen_fixed_fraction`.
Otherwise, I would suggest to revise the fractions you are using for grid
refinement.
Marc
On Wednesday, January 27, 2021 at 12:40:37 PM UTC-7 Karthi wrote:
> Hi Marc,
>
>
> Sorry for my delayed response, I was away for a couple of days. Thank you,
> your suggestions were very helpful and it worked.
>
> I have a couple of follow up questions in regard to mesh refinement.
>
>
> Question (1)
>
>
> I employ FE_Q(degree, 2) for all my sub-problems (and the same
> dof_handler).
>
>
> Hence, I have a std::vector, I created a container of
> estimated_error_per_cell
> for each solution and summed it all up in sum_estimater_error_per_cell,
> which is then used in the Kelly Error Estimator. I used fe.componet_mask to
> account for only the second component of my solution while calculating the
> error estimate.
>
>
> My objective is calculate a gradient error estimator as follow, where N
> represents solution.size() (i.e the number of sub-problems)
>
>
> \sum_{i}^{N} |\nabla \alpha_i|
>
>
>
> I don't know if the below code achieves this?
>
>
> std::vector> estimated_error_per_cell(num_index,Vector >(triangulation.n_active_cells()));
>
>
> Vector sum_estimated_error_per_cell(triangulation.n_active_cells
> ());
>
> FEValuesExtractors::Scalar alpha(1);
>
>
> for(unsigned int i=0; i < num_index; i++)
>
>KellyErrorEstimator::estimate(dof_handler,
>
> QGauss(fe.degree + 1),
>
> {},
>
> solution[i],
>
> estimated_error_per_cell[i],
>
> fe.component_mask(alpha));
>
>
> for(unsigned int i=0; i < num_index; i++)
>
> sum_estimated_error_per_cell += estimated_error_per_cell[i];
>
>
>
> GridRefinement::refine_and_coarsen_fixed_number(triangulation,
>
>
> sum_estimated_error_per_cell,
>
> 0.60,
>
> 0.40);
>
>
> Question (2)
>
>
> Please see the attachment of mesh refinement plots of the transient
> solution. As you can see, the initial mesh refinement at time zero makes
> sense but as the solution decays I was hoping the mesh would coarsen (but
> it refines further). I am clearly doing something wrong. I need some help
> in fixing this issue.
>
>
> Thank you!
>
>
> Karthi.
>
> On Mon, Jan 25, 2021 at 12:19 AM Marc Fehling wrote:
>
>> Hi Karthi,
>>
>> if you work on the same DoFHandler, one SolutionTransfer object is
>> sufficient.
>>
>> There are already member functions that take a container of solutions
>> like yours as a parameter. Have a look at
>> SolutionTransfer::prepare_for_coarsening_and_refinement
>>
>>
>> and SolutionTransfer::interpolate
>>
>> .
>>
>> Best,
>> Marc
>>
>> On Saturday, January 23, 2021 at 6:30:34 AM UTC-7 Karthi wrote:
>>
>>> Dear All,
>>>
>>> I have a fourth order parabolic equation, which is split into two second
>>> order equations. Hence I solve using components by
>>> declaring FESystem fe(FE_Q(degree,2). I have in total three such
>>> sub-problems, which are coupled to each other in a semi-implicit manner.
>>> Therefore I have a std::vector of solution for the entire system.
>>>
>>> std::vector> solution.
>>>
>>> In addition, I am employing mesh adaptivity. After estimating the error
>>> using Kelly, I would like to perform a solution transfer from old to new
>>> mesh. Do I need to create a std::vector of SolutionTransfer objects; one
>>> for each solution?
>>>
>>> The below code copied from step-26 seems to work. Is this the correct
>>> approach?
>>>
>>> std::vector> solution_trans(3,dof_handler);
>>>
>>> std::vector> previous_solution(num_index);
>>>
>>> for(unsigned int i=0; i < num_index; i++){previous_solution[i] =
>>> solution[i];}
>>>
>>> triangulation.prepare_coarsening_and_refine
