Re: [deal.II] Dealing with conflicting constraints without using VectorTools::interpolate_boundary_values()
I assign such a boundary condition at points by giving boundary values
to std::map boundary_values directly, and it looks
working well.
Yes, this is the correct approach.
However, for handling conflict between hanging node and
boundary condition with ConstraintMatrix, just doing something like:
typenameDoFHandler::active_cell_iterator cell = dof_handler.begin_active(),
endc =
dof_handler.end();
for(; cell!=endc; ++cell)
{
for(unsignedintvert = 0; vert < GeometryInfo::vertices_per_cell;
vert++)
{
if("certain condition")
{
constraints_update.add_line(cell2->vertex_dof_index(vert,I));
}
}
}
(Here constraints_update is an object of ConstraintMatrix).
does not seem working because the solver eventually gives the error:
[...]
Well, that means that your matrix is probably not invertible (although the
residual does get pretty small -- are you sure you don't just need more
iterations?). This *may* be because of your constraints, but it may also be
because of any number of other reasons. Have you tried this on a small
problem, say a 2x2 or 3x3 mesh, to make sure it really is because of the
constraints? For small problems, you can often figure out on a piece of paper
what the matrix should be, and compare with what you get from your program.
Best
W.
--
Wolfgang Bangerth email: bange...@colostate.edu
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 dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
Re: [deal.II] How to make G vector in lecture 21.65
Kyusik,
I'm studying Inhomogenous Dirichlet boundary conditions watching Lecture 21.65.
To practice this boundary condition, I'm trying to reproduce the result of
step-4(in which Inhomogenous Dirichlet boundary condition is used) in step-6.
As I learned in the lecture...
First, I assembled A and F with zero boundary condition by
template
void Step6::assemble_system ()
{
const QGauss quadrature_formula(3);
const RightHandSide right_hand_side;
FEValues fe_values (fe, quadrature_formula,
update_values| update_gradients |
update_quadrature_points | update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points= quadrature_formula.size();
FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell);
Vector cell_rhs (dofs_per_cell);
std::vector local_dof_indices (dofs_per_cell);
typename DoFHandler::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
cell_matrix = 0;
cell_rhs = 0;
fe_values.reinit (cell);
for (unsigned int q_index=0; q_indexget_dof_indices (local_dof_indices);
constraints.distribute_local_to_global (cell_matrix,
cell_rhs,
local_dof_indices,
system_matrix,
system_rhs);
}
}
What is the content of 'constraints' in your case?
And then, I made G_tilt by...
template
void Step6::assemble_G ()
{
const QGauss face_quadrature_formula(3);
FEFaceValues fe_face_values (fe, face_quadrature_formula,
update_values | update_gradients |
update_quadrature_points |
update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_face_q_points = face_quadrature_formula.size();
Vector local_rhs (dofs_per_cell);
G_tilt.reinit (dof_handler.n_dofs());
std::vector local_dof_indices (dofs_per_cell);
typename DoFHandler::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
local_rhs = 0;
for (unsigned int face_no=0; face_no::faces_per_cell;
++face_no)
if (cell->at_boundary(face_no))
{
fe_face_values.reinit (cell, face_no);
for (unsigned int q=0; q &p=fe_face_values.quadrature_point(q);
for (unsigned int i=0; iget_dof_indices (local_dof_indices);
for (unsigned int i=0; i
This is not correct. You are computing the G vector as
G_i = \int_{\partial\Omega} \varphi_i(x) g(x) ds
i.e., as a boundary integral. But G should be the vector that *interpolates*
g(x) on the boundary.
And I solved AU_0=F-AG_tilt by...
Vector tmp;
tmp.reinit(solution.size());
system_matrix.vmult(tmp,G_tilt);
system_rhs -= tmp;
solve ();
solution += G_tilt;
This looks correct.
Best
W.
--
Wolfgang Bangerth email: bange...@colostate.edu
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 dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
Re: [deal.II] Re: AMR , how to pass solution vector to refined mesh
On 03/17/2017 08:03 AM, Jaekwang Kim wrote:
After the first cycle, I am doing
template
voidnonlinear::refine_grid ()
{
Vector estimated_error_per_cell (triangulation.n_active_cells());
KellyErrorEstimator::estimate (dof_handler,
QGauss(3),
typenameFunctionMap::type(),
solution,
estimated_error_per_cell);
GridRefinement::refine_and_coarsen_fixed_number (triangulation,
estimated_error_per_cell,
0.3, 0.0);
triangulation.prepare_coarsening_and_refinement ();
SolutionTransfer solution_transfer(dof_handler);
solution_transfer.prepare_for_coarsening_and_refinement(solution);
triangulation.execute_coarsening_and_refinement(); // I think n_dof() has
increased here. You need to check this if necessary
dof_handler.distribute_dofs(fe);
previous_solution.reinit (dof_handler.n_dofs());
solution_transfer.interpolate(solution, previous_solution);
DynamicSparsityPattern dsp(dof_handler.n_dofs());
DoFTools::make_sparsity_pattern (dof_handler, dsp);
sparsity_pattern.copy_from(dsp);
solution.reinit (dof_handler.n_dofs());
system_matrix.reinit (sparsity_pattern);
system_rhs.reinit (dof_handler.n_dofs());
}
I am sure that I am reinitialize matrixes...
but I am not sure about whether the redline is doing enough for
reinitialization of sparsity pattern
It looks like it, but it seems like it's not working, and just seeing the
little code snippet is not enough for me to tell for sure. You'll have to
debug what is going on. My first stab would be: At the beginning of your
assembly function, put something like
Assert (system_matrix.n() == dof_handler.n_dofs(), ExcInternalError());
to verify that indeed your matrix has the right size.
Best
W.
--
Wolfgang Bangerth email: bange...@colostate.edu
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 dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
Re: [deal.II] Transfer local cell data to global cell vector
I need to do some post-processing on computed solution data (say Vector Solution(dof_handler.n_dofs()) is the solution vector). For post-processing field (say Vector Curvature(triangulation.n_dofs()) is the scalar field to be computed) I decided to have one entry per cell. Once I calculate this entry in 'cell iterator loop' I have to transfer it to global vector Curvature based on some index of the cell. I am just drawing the analogy from what we do during the assembling of system_matrix and/or system_rhs. There we use cell iterator functionality get_dof_indices() to get the global indices of DoFs that are of on the current cell. But I am not able to find analogous method to get an cell-wise index (we can't use cell iterator return value since it is not unsigned int) which can be used to transfer cell data to the global vector Curvature(triangulation.n_dofs()). Please help me to figure it out. You can use cell->active_cell_index(). Secondly, I want to use this global vector Curvature(triangulation.n_dofs()) in a DataOut object essentially to add this vector to the object via add_data_vector(Curvature, "curvature") method. I believe (read in another question thread in the group) the DataOut object will not have any problems accepting a vector which has size triangulation.n_dofs() but not dof_handler.n_dofs() since I initialize the DataOut object with dof_handler? I meant that it would judiciously distinguish this vector from other vectors which have size dof_handler.n_dofs() and output it correctly? It can take both kinds of vectors. In your case, rather than converting the cell-based vector into a DoF-based vector, stick with the former and let DataOut output this. Best W. -- Wolfgang Bangerth email: bange...@colostate.edu 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 dealii+unsubscr...@googlegroups.com. For more options, visit https://groups.google.com/d/optout.
Re: [deal.II] How to make G vector in lecture 21.65
Dr. Bangerth,
Thank you very much for your kind reply.
> What is the content of 'constraints' in your case?
>
> DoFTools::make_hanging_node_constraints (dof_handler,
constraints);
VectorTools::interpolate_boundary_values (dof_handler,
0,
ZeroFunction(),
constraints);
constraints.close ();
DynamicSparsityPattern dsp(dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler,
dsp,
constraints,
/*keep_constrained_dofs = */ false);
sparsity_pattern.copy_from(dsp);
system_matrix.reinit (sparsity_pattern);
}
As you said in the lecture, I apply the ZeroFunction() for
A*U_0=F-A*G_tilt
>
> This is not correct. You are computing the G vector as
>G_i = \int_{\partial\Omega} \varphi_i(x) g(x) ds
> i.e., as a boundary integral. But G should be the vector that
> *interpolates*
> g(x) on the boundary.
>
> template
void Step6::assemble_G ()
{
const QGauss face_quadrature_formula(3);
FEFaceValues fe_face_values (fe, face_quadrature_formula,
update_values | update_gradients |
update_quadrature_points |
update_JxW_values);
const unsigned int dofs_per_face = fe.dofs_per_face;
const unsigned int n_face_q_points = face_quadrature_formula.size();
Vector local_rhs (dofs_per_face);
G_tilt.reinit (dof_handler.n_dofs());
std::vector local_dof_indices (dofs_per_face);
typename DoFHandler::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
local_rhs = 0;
for (unsigned int face_no=0;
face_no::faces_per_cell;
++face_no)
if (cell->at_boundary(face_no))
{
fe_face_values.reinit (cell, face_no);
for (unsigned int q=0; q &p=fe_face_values.quadrature_point(q);
for (unsigned int i=0; iface(face_no)->get_dof_indices (local_dof_indices);
}
for (unsigned int i=0; i const
unsigned int dofs_per_face = fe.dofs_per_face;
2. Vector local_rhs(dofs_per_cell); -> Vector
local_rhs (dofs_per_face);
3.std::vector local_dof_indices (dofs_per_cell);
-> std::vector local_dof_indices (dofs_per_face);
4.
for (unsigned int i=0; i
for (unsigned int i=0; iget_dof_indices (local_dof_indices);
for (unsigned int i=0; i
cell->face(face_no)->get_dof_indices (local_dof_indices);
for (unsigned int i=0; ihttp://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 dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
[deal.II] Re: Transfer local cell data to global cell vector
Thankyou! that works :) On Saturday, 18 March 2017 12:03:55 UTC+5:30, Sumedh Yadav wrote: > > Hello, > I need to do some post-processing on computed solution data (say > Vector Solution(dof_handler.n_dofs()) is the solution vector). For > post-processing field (say Vector Curvature(triangulation.n_dofs()) > is the scalar field to be computed) I decided to have one entry per cell. > Once I calculate this entry in 'cell iterator loop' I have to transfer it > to global vector Curvature based on some index of the cell. I am just > drawing the analogy from what we do during the assembling of system_matrix > and/or system_rhs. There we use cell iterator functionality > get_dof_indices() to get the global indices of DoFs that are of on the > current cell. But I am not able to find analogous method to get an > cell-wise index (we can't use cell iterator return value since it is not > unsigned int) which can be used to transfer cell data to the global vector > Curvature(triangulation.n_dofs()). Please help me to figure it out. > Secondly, I want to use this global vector > Curvature(triangulation.n_dofs()) in a DataOut object essentially to add > this vector to the object via add_data_vector(Curvature, "curvature") > method. I believe (read in another question thread in the group) the > DataOut object will not have any problems accepting a vector which has size > triangulation.n_dofs() but not dof_handler.n_dofs() since I initialize the > DataOut object with dof_handler? I meant that it would judiciously > distinguish this vector from other vectors which have size > dof_handler.n_dofs() and output it correctly? > -- 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 dealii+unsubscr...@googlegroups.com. For more options, visit https://groups.google.com/d/optout.
