Re: [deal.II] Computing the curl of a solution vector field obtained from Nedelec elements
Sebastian, > I am currently experimenting with curl-conforming elements and their > application in electromagnetism. I am using the Nedelec element FE_Nedelec > and > I have difficulties to understand the output and some other features: > > 1. I realized that VectorTools::project only works if QGauss(p) with p > >= 2 is used. AssertNoZerosOnDiagonal fails in precondition_SSOR. If I am > using an FE_System(FE_Q(1), dim) as the base element, everything > is fine. I don't know if this qualifies as a bug. But there is no > indication in the documentation of certain requirements on the quadrature > applied. There is a statement, but you're right that it isn't quite clear. I've proposed this here: https://github.com/dealii/dealii/pull/7587 > 2. I tried to output the curl of the vector, A, to a vtk file and used a > class derived from DataPostprocessor as suggested above. The results > do not look very promising if I compare them with the output of an > analytic curl of A, see screenshots attached. However the output of the > potential itself works as expected. It's hard to see anything in these figures if you don't know what it represents :-) I looked at the code briefly and it looks correct (though I'm not sure about the sign of each entry of the curl.) Try this on a simple case first where you take a simple vector-valued function (say, one whose curl is a constant vector), project it onto the mesh, and then output both the vector potential and the curl. > 3. Because of the problem related to the output of the curl I tried to > compute some errors. The projection error of the function value itself is > converging pretty well. I have tried to compute the L2-error of curl(A_h) > - curl(A) where A_h is the discrete approximation of A and A is the exact > potential. The curl error is not converging if I am refining the mesh. But > I am not sure if there is a mistake in my computation of the > cellwise_error. Does the potential A_h itself converge at the correct order? If yes, then its curl should converge at one order less. But since it's a derived quantity, make sure first that the quantity from which it is derived (i.e., the primary variable A_h) converges correctly and with the correct order. If it does and the curl doesn't, then you know that the problem lies with the way you compute the curl, not the way you compute A_h. 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] Usage of FEFieldFunction.vector_value_list on a parallel::distributed::Triangulation
On 1/8/19 12:43 AM, 'Maxi Miller' via deal.II User Group wrote: > > I assume it tries to find the point, which is not available due to being > within a cell located in another thread. I assume I have to modify the cell > iterator in that case. Is that correct? You mean owned by some other MPI process? I don't know whether you interpretation is correct. The exception is found here: // The first time around, we check for vertices in the hint_cell. If // that does not work, we set the cell iterator to an invalid one, and // look for a global vertex close to the point. If that does not work, // we are in trouble, and just throw an exception. // // If we got here, then we did not find the point. If the // current_cell.state() here is not IteratorState::valid, it means that // the user did not provide a hint_cell, and at the beginning of the // while loop we performed an actual global search on the mesh // vertices. Not finding the point then means the point is outside the // domain. AssertThrow(current_cell.state() == IteratorState::valid, ExcPointNotFound(p)); Do you think you can come up with a small testcase? 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: Implement of adaptive mesh refinement cannot work in time-dependent problem.
Chucui,
> And I think what you means is the same as it above beforehand, but maybe not
> now. The dsp (BlockDynamicSparsityPattern) will be destroyed and cannot exist
> in all lifetime of BlockSparseMatrix, so we define : BlockSparsityPattern
> sparsity_pattern at the top, (this is what step-22 does) and let it copy
> dsp
> But now, we know that BlockSparsityPattern sparsity_pattern also will be
> destroyed and cannot exist in all lifetime of BlockSparseMatrix. So I don't
> know how to correct it.
> And I cannot implement what you said "You'll first have to break that
> dependence before you destroy the sparsity pattern".
What I'm trying to say is this: Let's say you have code like this:
{
SparseMatrixm;
SparsityPattern sp;
...fill sp, for example by copying from a DynamicSparsityPattern...
m.reinit (sp); // m now points to sp
...some more code...
}
At the '}' brace, the compiler destroys objects in the reverse order as they
were created. So it first destroys 'sp'. But 'm' still points to it and so you
get an error: since 'sp' is still needed by the variable 'm', you can't
destroy 'sp' yet. If you want to destroy 'sp', you first have to tell 'm' to
not use 'sp' any more, for example by using code such as this:
{
SparseMatrixm;
SparsityPattern sp;
...fill sp, for example by copying from a DynamicSparsityPattern...
m.reinit (sp); // m now points to sp
...some more code...
m.reinit(); // ***
}
In the marked line, you tell 'm' to not use 'sp' any more. So at the '}'
brace, 'sp' can be destroyed now because no other variable is using it any
more and you will not get the error any more. Next 'm' is destroyed, but that
is no longer of concern to us.
> 2. As your suggestions, I change my code :
> |
> template
> void
> StokesProblem::refine_mesh (const unsigned int max_grid_level,
> const unsigned int min_grid_level)
> {
> Vector estimated_error_per_cell
> (triangulation_active.n_active_cells());
> FEValuesExtractors::Scalar phase(0);
> KellyErrorEstimator::estimate (dof_handler,
> QGauss(degree+1),
> typename FunctionMap::type(),
> old_solution,
> estimated_error_per_cell,
> fe.component_mask(phase));
> GridRefinement::refine_and_coarsen_fixed_number (triangulation_active,
>
> estimated_error_per_cell,
> 0.2,0.1);
> std::cout << "triangulation.n_levels" << triangulation_active.n_levels() <<
> std::endl;
> if (triangulation_active.n_levels() > max_grid_level)
> for (typename Triangulation::active_cell_iterator
> cell = triangulation.begin_active(max_grid_level);
> cell != triangulation.end(); ++cell)
> cell->clear_refine_flag ();
> std::cout << "clear-refine" << triangulation_active.n_levels() << std::endl;
> if (triangulation_active.n_levels() < min_grid_level)
> for (typename Triangulation::active_cell_iterator
> cell = triangulation.end_active(min_grid_level);
> cell != triangulation.end(); ++cell)
> cell->clear_coarsen_flag ();
> std::cout << "clear-coarsen" << triangulation_active.n_levels() << std::endl;
> triangulation_active.execute_coarsening_and_refinement ();
> std::cout << "c-and-r" << triangulation_active.n_levels() << std::endl;
> }
> |
>
> Maybe I need to chage all"triangulation" in my originial
> StokesProblem::refine_mesh (const unsigned int max_grid_level, const
> unsigned int min_grid_level) into "triangulation_active", this function can
> run(but i don't know if it is right), and for implement periodic boundary
> condition, I use the "triangulation" and "triangulation_active" :
I don't know what the difference between triangulation and
triangulation_active is. It's really hard to debug a code you don't have or
know :-) But I will point out that something like this here:
if (triangulation_active.n_levels() > max_grid_level)
for (typename Triangulation::active_cell_iterator
cell = triangulation.begin_active(max_grid_level);
looks suspicious: In the 'if' condition, you are referencing
triangulation_active, but then when you are asking for a concrete cell, you
are referencing triangulation.begin_active(max_grid_level). This doesn't look
right.
> I don't know what I should use in whole code is "triangulation" or
> "triangulation_active" ?
A question like this suggests to me that you're not quite clear about what
these two objects are and where they should be used. That sounds like you
don't quite know what your code does or *should* do. Do you have colleagues
with whom you could talk through the
[deal.II] deal.II Newsletter #62
Hello everyone! This is deal.II newsletter #62. It automatically reports recently merged features and discussions about the deal.II finite element library. ## Below you find a list of recently proposed or merged features: #7580: Use more flags in [[fallthrough]] checks (proposed by masterleinad) https://github.com/dealii/dealii/pull/7580 #7578: Compute Point Location new function (proposed by GivAlz) https://github.com/dealii/dealii/pull/7578 #7575: Use range-based for loops II (proposed by masterleinad) https://github.com/dealii/dealii/pull/7575 #7574: Check [[fallthrough]] instead of assuming it when C++17 support is enabled (proposed by masterleinad; merged) https://github.com/dealii/dealii/pull/7574 #7573: Require OpenMP for CUDA support (proposed by masterleinad; merged) https://github.com/dealii/dealii/pull/7573 #7572: Fixed indices in FESeries::Fourier::set_k_vectors(). (proposed by marcfehling) https://github.com/dealii/dealii/pull/7572 #7571: Replace size_t by std::size_t (proposed by masterleinad; merged) https://github.com/dealii/dealii/pull/7571 #7570: Use range-based for loops in examples and expand_instantiations (proposed by masterleinad; merged) https://github.com/dealii/dealii/pull/7570 #7569: Add test and support for non-mpi to GridTools::build_global_tree (proposed by GivAlz; merged) https://github.com/dealii/dealii/pull/7569 #7568: Some minor cleanups in GridTools. (proposed by bangerth; merged) https://github.com/dealii/dealii/pull/7568 #7567: Update some documentation. (proposed by bangerth; merged) https://github.com/dealii/dealii/pull/7567 #7565: Use correct data type for a random number. (proposed by bangerth; merged) https://github.com/dealii/dealii/pull/7565 #7564: Optimizing NonMatching::create_coupling_mass_matrix (proposed by GivAlz; merged) https://github.com/dealii/dealii/pull/7564 #7563: Optimizing NonMatching::create_coupling_sparsity_pattern (proposed by GivAlz; merged) https://github.com/dealii/dealii/pull/7563 #7562: Fix nvcc warnings (proposed by masterleinad; merged) https://github.com/dealii/dealii/pull/7562 #7474: Fix bugs in extract_boundary_mesh (proposed by starki0815; merged) https://github.com/dealii/dealii/pull/7474 #7440: Disallow using rvalue references for Subscriptor::(un)?subscribe (proposed by masterleinad; merged) https://github.com/dealii/dealii/pull/7440 #7303: Allow using CUDA-aware MPI (proposed by masterleinad; merged) https://github.com/dealii/dealii/pull/7303 #7020: Add HDF5 classes (proposed by dangars; merged) https://github.com/dealii/dealii/pull/7020 #7001: Adding global tree of bounding boxes (proposed by GivAlz; merged) https://github.com/dealii/dealii/pull/7001 ## And this is a list of recently opened or closed discussions: #7581: copy_triangulation removes limit_level_difference_at_vertices flag (opened) https://github.com/dealii/dealii/issues/7581 #7579: Error "lambda-expression in unevaluated context" in AssertDimension (opened and closed) https://github.com/dealii/dealii/issues/7579 #7577: Should GridTools::extract_used_vertices return also artificial vertices? (opened) https://github.com/dealii/dealii/issues/7577 #7576: Behaviour of find active cell around point with artificial cells (opened) https://github.com/dealii/dealii/issues/7576 #7566: Rename SolutionTransfer::prepare_serialization (opened) https://github.com/dealii/dealii/issues/7566 #7551: Possible extension of SmoothnessEstimator::legendre_coeff_decay (closed) https://github.com/dealii/dealii/issues/7551 #7467: Bug in GridGenerator::extract_boundary_mesh (closed) https://github.com/dealii/dealii/issues/7467 #7389: deal.II 9.0 not working on CentOS 6 (closed) https://github.com/dealii/dealii/issues/7389 #7352: Underlinkage in quick tests with GSL (closed) https://github.com/dealii/dealii/issues/7352 A list of all major changes since the last release can be found at https://www.dealii.org/developer/doxygen/deal.II/changes_after_8_5_0.html. Thanks for being part of the community! Let us know about questions, problems, bugs or just share your experience by writing to dealii@googlegroups.com, or by opening issues or pull requests at https://www.github.com/dealii/dealii. Additional information can be found at https://www.dealii.org/. -- 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.
[deal.II] Re: copy_triangulation removes limit_level_difference_at_vertices flag
Done: https://github.com/dealii/dealii/issues/7581 W dniu wtorek, 8 stycznia 2019 23:42:06 UTC+1 użytkownik Daniel Arndt napisał: > > Michal, > > not setting `limit_level_difference_at_vertices` when copying a > parallel::distributed::Triangulation can be considered an oversight/bug. > Do you want to create a PR for fixing this? > > Best, > Daniel > -- 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] Computing the curl of a solution vector field obtained from Nedelec elements
You shouldn't use FE_Nedelec with non-standard meshes like hyper_ball. Try FE_NedelecSZ instead. Read their docs for more information. Best Alexander -- 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.
