Dear Dr. Wolfgang Bangerth, Thanks for your response. I am quite interested in how to implement projection.
In fact, I have tried both methods but did not get the correct result using
projection:
- Interpolate the boundary data: It works quite well for adaptivity with
both 2 and 3 dimensions. Though expecting extra error will be introduced, I
observed that it is actually comparable to the case with projected boundary
data (I do the experiment through global refinement).
- Project the boundary data: I tried to project the boundary data by
commenting out the project_boundary_values and adding attached code in
step-51:
if (task_data.trace_reconstruct == false)
{
...
cell->get_dof_indices(task_data.dof_indices);
my code
}else{
...
}
where face_mass is the size of fe.dofs_per_cell x fe.dofs_per_cell matrix,
and both dirichlet_rhs and proj_dirichlet_rhs are fe.dofs_per_cell vector.
I just solve the linear system on the boundary faces, assign the projected
data to the cell_vector, and zero out the off-diagonal terms. In this
procedure, I assume constraints can handle hanging nodes properly. Whatever
values I assign to the corresponding row will be replaced by some algebraic
function that will be used to maintain the continuity via data from end
support points.
The code can run but the linear solver for solving trace unknowns requires
a significant iteration number (will explode if using global refinement
option). In addition, the L2 error norm is too high compared to the
original setting. Hence, I may do something wrong. Could you give me some
suggestions or comments on it?
Best Regards,
Jau-Uei Chen
Wolfgang Bangerth 在 2022年3月20日 星期日下午9:49:49 [UTC-5] 的信中寫道：
>
> > I am currently working on a HDG solver and use step-51 as a reference. I
> found
> > that the example does not work for a 3d setting along with adaptivity.
> It
> > seems that an error is from "VectorTools::project_boundary_values", the
> error
> > message I got is:
> >
> > An error occurred in line <742> of file
> >
> </org/groups/ceo/chenju/candi/install/tmp/unpack/deal.II-v9.2.0/include/deal.II/numerics/vector_tools_boundary.templates.h>
>
>
> > in function
> > void dealii::VectorTools::internal::do_project_boundary_values(const
> > M_or_MC<dim, spacedim>&, const DoFHandlerType<dim, spacedim>&, const
> > std::map<unsigned int, const dealii::Function<spacedim, number>*>&,
> const
> > Q_or_QC<(dim - 1)>&,std::map<unsigned int, number>&,
> std::vector<unsigned
> > int>) [with int dim = 3; int spacedim = 3; DoFHandlerType =
> > dealii::DoFHandler; M_or_MC = dealii::Mapping; Q_or_QC =
> dealii::Quadrature;
> > number = double]
> > The violated condition was:
> > level == cell->level()
> > Additional information:
> > The mesh you use in projecting boundary values has hanging nodes at
> the
> > boundary. This would require dealing with hanging node constraints when
> > solving the linear system on the boundary, but this is not currently
> implemented.
> >
> > The version of dealii I use is 9.2.0. I am wondering if there is any way
> to
> > resolve this error.
>
> Not easily. Basically, what is happening is that for *projection* of
> boundary
> values, you need to build and solve a linear system of only those degrees
> of
> freedom that are located on the boundary. The way this is currently done
> is
> unable to deal with hanging node constraints.
>
> There are approaches to deal with this, but they are (as far as I know)
> not
> currently implemented and so even newer versions of deal.II will not
> support
> this. There are two options for you:
> * You can use *interpolation* (rather than projection) of boundary values
> * You can implement the missing functionality. We'd be happy to guide you
> on
> how this can be done, but it will involve a certain amount of work.
>
>
> > Besides this error, I also have a hard time investigating
> > the vertice initiated by FE_FaceQ and would also like to seek some
> > suggestions. Here is my code modified from step-51:
> >
> > template <int dim>
> > void HDG<dim>::assemble_system_one_cell(
> > const typename DoFHandler<dim>::active_cell_iterator &cell,
> > ScratchData & scratch,
> > PerTaskData & task_data)
> > {
> > for( unsigned int face = 0; face < GeometryInfo< dim
> > >::faces_per_cell; ++face ){
> > if(cell->face(face)->at_boundary() &&
> > (cell->face(face)->boundary_id() == 0)){
> > std::cout << "node points" << std::endl;
> > for (const auto v: cell->face(face)->vertex_indices()){
>
> This uses a function that was introduced only in deal.II 9.3, I believe.
> But
> you can equivalently implement this line using
> for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
>
>
> 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].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/42780ee0-4a11-4f51-a00d-0093d27d9dd7n%40googlegroups.com.
project boundary data.pdf
Description: Adobe PDF document
