Re: [deal.II] Using solution on from one FE problem as a boundary condition for another

[Daniel Arndt]

Ernesto,

Imposing strong boundary conditions should not be overly difficult in this
case.
Just as Wolgang said you should have a look at
interpolate_boundary_values() and replace the Function object evaluations
on each cell by the function values given through
a FEValues object initialized for your first solution. Assuming that you
are using interpolatory ansatz functions, this interpolation has to happen
at the support points of your target finite element.
In particular, these support points should be the quadrature points for the
FEValues object.

Best,
Daniel

Am Di., 14. Jan. 2020 um 01:50 Uhr schrieb Ernesto Ismail <
ernesto.ism...@gmail.com>:

> Dear Prof. Bangerth,
>
> Thank you for your reply.
>
> I would ideally have liked to impose the boundary conditions strongly but
> it looks like doing so will not be straight-forward. I will do some
> investigations into applying the boundary conditions weakly and attempt
> that first.
>
> One approach I had tried was to use Lobatto quadrature and use a
> quadrature point history data structure to transfer the solution at the
> boundary. The problem with this approach was that the differing element
> types lead to some strange behaviours along the boundary.
>
> Thanks again,
> Ernesto
>
> On Thursday, 9 January 2020 23:43:27 UTC+2, Wolfgang Bangerth wrote:
>>
>>
>> Ernesto,
>>
>> > I am currently working with a staggered coupled system where I am
>> solving
>> > several finite element problems in sequence. Up until now my boundary
>> > conditions have been fairly straightforward and I've been able to get
>> > meaningful results from my work.
>> >
>> > The finite element schemes I use for each step are fairly different,
>> with some
>> > being discontinuous and of varying orders.
>> >
>> > I am now in a position where I need to use the solution of one of my
>> steps as
>> > the boundary condition for the next (I actually need to combine two of
>> them
>> > with a simple arithmetic function, but baby steps for now). I presume I
>> need
>> > to use some sort of projection to achieve this, but I'm not exactly
>> sure how
>> > to go about it. I thought to try and use the project_boundary_values()
>> > function to do what I need, but I am unsure how I would structure
>> > the boundary_functions argument.
>>
>> Do you need to enforce these boundary conditions strongly? If you can
>> enforce
>> them weakly, then the solution of one problem only enters the weak form
>> of
>> another problem (in that case through a boundary integral) and that is
>> really
>> not very different from the way we do this in other problems where we
>> have
>> multiple DoFHandlers -- e.g., in step-31. (Though my usual advice still
>> holds
>> true that I think everything becomes simpler if you pack everything into
>> one
>> DoFHandler).
>>
>> It gets more complicated if you want to impose strong boundary
>> conditions.
>> You've already found the FEFieldFunction function, but it's complicated
>> and
>> slow. A simpler way would be if you looked at the implementation of the
>> interpolate_boundary_values() function and identified where it is asking
>> the
>> Function object for its values. That's the place where you'd need to
>> somehow
>> combine your existing solutions into something else and turn that into a
>> constraint.
>>
>> Best
>>   W.
>>
>> --
>> 
>> Wolfgang Bangerth  email: bang...@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.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/832bb83a-7e04-4243-9524-64db1bfc9345%40googlegroups.com
> 
> .
>

-- 
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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/CAOYDWbKY2%3D%2BBdWX4mh7guYEROFPT7mX8QR%3D4ED4KWymi_cD7Aw%40mail.gmail.com.

Re: [deal.II] discontinous contour over elements

[David Eaton]

Dear Bruno,

Thank you for your suggestions. I am going to take a look at Lethe and compare 
with my implementation. In stabilized formulation, I used quadrilateral 
element, instead of P2 P1 Taylor-Hood element. The used element is only C0 
element. I also did not expect such a discontinuity between elements. Although 
I use P2 P1 Taylor-Hood element without stabilization terms, the discontinuity 
is still there. Probably I made mistakes somewhere  in setup. I also suspect 
that my solution is not converged. After taking a small relative tolerance 
1e-8, the discontinuity still appears. As Professor Wolfgang suggested, I am 
currently checking the convergence rate of this formulation. Thank you for your 
suggestions. If I cannot resolve this issue, I will update in the forum again.

Best
D.

From: dealii@googlegroups.com  on behalf of Bruno 
Blais 
Sent: Tuesday, January 14, 2020 9:53 PM
To: deal.II User Group 
Subject: Re: [deal.II] discontinous contour over elements

Dear David,
How are you calculating the vorticity?
As Wolfgang and Praveen have mentioned, if you are using the DataPostProcessor, 
then this will use your shape functions to calculate the vorticity. However, 
your P2-P1 elements are only C0 continuous. Consequently, your vorticity can 
possibly be inherently discontinuous at the element edges. However, I am 
surprised that you obtain such strong discontinuity. In our code based on 
deal.ii (https://github.com/lethe-cfd/lethe) we implement stabilized 
formulations for the NS equations and the vorticity results for such cases are 
very smooth (even when represented using discontinuous shape functions.

Have you established the convergence of your code using manufactured solution? 
This is the first thing I would suggest. You can look at the applications_tests 
of lethe for some examples of easy manufactured solutions for the 
Incompressible Navier-Stokes equations. There are also common books that treat 
this issue (for instance : 
https://www.amazon.ca/Verification-Validation-Scientific-Computing-Oberkampf/dp/0521113601)

Please feel free to reach out if you have more questions.
Best
Bruno




On Monday, 13 January 2020 20:29:53 UTC-5, David Eaton wrote:
Thanks Wolfgang and Praveen for providing suggestions. I have tried to 
debugging this code for a while. I have attached this simple code on this 
email. I followed the instructions in tutorial closely. Hopefully, anyone could 
give some suggestions.

Best
D.

From: dea...@googlegroups.com  on behalf of Wolfgang 
Bangerth 
Sent: Tuesday, January 14, 2020 6:24 AM
To: dea...@googlegroups.com 
Subject: Re: [deal.II] discontinous contour over elements

On 1/12/20 9:17 PM, David Eaton wrote:
> My inflow condition is uniform. This formulation and mesh is tested in a
> simple C++ code without library. The  large mesh near the inflow does not give
> this problem.
> Yes. I am using C0 element. I did calculation using tecplot. However, the
> results from a my C++ code does not give this issue either. Just now, I check
> the formulation again. Although I use Q2Q1 Taylor-Hood element without any
> stabilization, these issues are still happening.

David -- we don't really know what formulation you are using, how you are
implementing it, what you are comparing against, and a number of other factors.

If you have a formulation that computes u,p, and you are plotting the
vorticity, you need to expect that the isocontours are discontinuous for the
reasons Praveen already stated. If you are getting results that make no sense
to you, then the first step would be to ensure that your program is converging
as expected. To do this, choose a solution that you know and compute the error
norm; then ensure that the program yields error norms that decrease as
expected with mesh refinement.

Best
  W.

--

Wolfgang Bangerth  email: bang...@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 dea...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/e1b55e95-a139-d62e-e786-61393cf84dc7%40colostate.edu.

--
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 

Re: [deal.II] Re: KDTree question

[Wolfgang Bangerth]

> |
> |template
> classA
> {
> public:
>      A(std::vector>);
> constKDTree_kdtree()const;
> doubledoing_something();
> 
> private:
> KDTreekdtree;
> };
> |
> |
> 
> and in my implementation i wrote
> |
> template
> A::A(std::vector>)
> {
>      kdtree.set_points(points);
> }
> 
> template
> constKDTree&
> A::get_kdtree()const
> {
> returnkdtree;
> }
> 
> template
> double
> A::doing_something()
> {
> constauto&_v_kdtree =this->get_kdtree();
> 
> }
> 
> |

Thus looks right to me.


> But everytime i called the get_kdtree function, i got the segmentation fault( 
> you may see the snippet in doing_something() function. Any idea to resolve 
> this problem?

We just had another thread on the mailing list about segmentation faults and 
how one debugs them. I don't think that your get_kdtree() function is wrong. 
Something must be wrong in how you set up the object that contains it. You 
will need to find out which pointer it is that creates the segmentation fault; 
if you know this part, it is typically easy to understand what is happening.

My suspicion is that in your case, the 'this' pointer is a nullptr.

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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/ae158c1f-a025-8d97-c5b0-fd8f7365d070%40colostate.edu.

[deal.II] deal.II Newsletter #106

[Rene Gassmoeller]

Hello everyone!

This is deal.II newsletter #106.
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:

#9318: converts l_2 to LaTeX code (proposed by krishnakumarg1984) 
https://github.com/dealii/dealii/pull/9318

#9315: In step 6, removes the documentation sentences relating to quadrature 
formula degree since already done in Step-5 (proposed by krishnakumarg1984; 
merged) https://github.com/dealii/dealii/pull/9315

#9313: In step6, removes outdated sentence on destructor (proposed by 
krishnakumarg1984; merged) https://github.com/dealii/dealii/pull/9313

#9311: Add a few links to step-6. (proposed by bangerth; merged) 
https://github.com/dealii/dealii/pull/9311

#9310: Disentangle the dependency graph some more. (proposed by bangerth; 
merged) https://github.com/dealii/dealii/pull/9310

#9309: FESeries: Added features to precalculate all matrices and serialize 
them. (proposed by marcfehling) https://github.com/dealii/dealii/pull/9309

#9304: In step 6 doc, fixes a really minor grammar/typo (proposed by 
krishnakumarg1984; merged) https://github.com/dealii/dealii/pull/9304

#9303: In the doc of the point() function, replaces monospace l_2 with its aTeX 
equivalent. Fixes #9302 (proposed by krishnakumarg1984) 
https://github.com/dealii/dealii/pull/9303

#9301: GridGenerator::half_hyper_shell warning (proposed by tjhei) 
https://github.com/dealii/dealii/pull/9301

#9300: In step 6 results doc, move the header file inclusion sentence to 
reflect that we need it for Sparse_ILU and not Jacobi preconditioner (proposed 
by krishnakumarg1984) https://github.com/dealii/dealii/pull/9300

#9298: Clarify GridGenerator::hyper_shell in 3D (proposed by masterleinad) 
https://github.com/dealii/dealii/pull/9298

#9297: Make GridTools::delete_duplicate_vertices() faster. (proposed by 
drwells) https://github.com/dealii/dealii/pull/9297

#9295: fixes stale documentation about boundary_object for hyper_ball function 
(proposed by krishnakumarg1984; merged) 
https://github.com/dealii/dealii/pull/9295

#9293: minor grammar fixes to step-6 tutorial doc (proposed by 
krishnakumarg1984; merged) https://github.com/dealii/dealii/pull/9293

#9291: In the doc for KellyErrorEstimator class, trying to fix couple of 
doxygen links (proposed by krishnakumarg1984; merged) 
https://github.com/dealii/dealii/pull/9291

#9290: fix write_vtu_in_parallel invalid file if partly empty (proposed by 
tjhei; merged) https://github.com/dealii/dealii/pull/9290

#9289: Remove filtered_iterator include in step-18 (proposed by peterrum; 
merged) https://github.com/dealii/dealii/pull/9289

#9288: Update constraints documentation module to reflect recent updates to 
step-6 (proposed by krishnakumarg1984) 
https://github.com/dealii/dealii/pull/9288

#9287: Replace the Dataout::first_cell/next_cell() mechanism (proposed by 
bangerth) https://github.com/dealii/dealii/pull/9287

#9286: Introduce \jump and \average commands to be used in the documentation. 
(proposed by bangerth; merged) https://github.com/dealii/dealii/pull/9286

#9285: Avoid gcc-4.9.4 error (proposed by masterleinad; merged) 
https://github.com/dealii/dealii/pull/9285

#9284: Fix \coloneq again. (proposed by bangerth; merged) 
https://github.com/dealii/dealii/pull/9284

#9283: Remove tensor.cc (proposed by masterleinad; merged) 
https://github.com/dealii/dealii/pull/9283

#9282: update .gitattributes (proposed by tjhei; merged) 
https://github.com/dealii/dealii/pull/9282

#9281: use github CI for OSX (proposed by tjhei; merged) 
https://github.com/dealii/dealii/pull/9281

#9280: Multigrid Transfer cleanup (proposed by tjhei; merged) 
https://github.com/dealii/dealii/pull/9280

#9278: implement SingleLevelDataOut (proposed by tjhei) 
https://github.com/dealii/dealii/pull/9278

#9277: Fix the dealcoloneq command. (proposed by bangerth; merged) 
https://github.com/dealii/dealii/pull/9277

#9276: Let latex figure out spacing itself. (proposed by bangerth; merged) 
https://github.com/dealii/dealii/pull/9276

#9275: replace the use of 'he' with 'it' to refer to neighbouring cell in 
tria.cc (proposed by krishnakumarg1984; merged) 
https://github.com/dealii/dealii/pull/9275

#9274: fixes a keyboard fat-fingering typo wherever 'he' in incorrectly used 
instead of 'the' (proposed by krishnakumarg1984; merged) 
https://github.com/dealii/dealii/pull/9274

#9273: Rewrite wording in headers/docs such that we now have gender neutral 
wording throughout (proposed by krishnakumarg1984; merged) 
https://github.com/dealii/dealii/pull/9273

#9272: Make DoFTools::count_dofs_per_{component,block} return its data. 
(proposed by bangerth) https://github.com/dealii/dealii/pull/9272

#9270: Step-3 results/extension documentation improvements and fixes. Closes 
#9227 (proposed by krishnakumarg1984; merged) 
https://github.com/dealii/dealii/pull/9270

#9268: Use AssertIndexRange some more 

[deal.II] Re: KDTree question

[A.Z Ihsan]

to be precise
template
class A
{
public: 
A(std::vector> );
const KDTree & get_kdtree() const;
double doing_something();

private:
KDTree kdtree;
};

and in my implementation i wrote
template 
A::A(std::vector> )
{
kdtree.set_points(points);
}

template 
const KDTree & 
A::get_kdtree() const
{
return kdtree;
}

template 
double
A::doing_something()
{
const auto& _v_kdtree = this->get_kdtree();

}


But everytime i called the get_kdtree function, i got the segmentation 
fault( you may see the snippet in doing_something() function. Any idea to 
resolve this problem?



On Tuesday, January 14, 2020 at 1:08:32 PM UTC+1, A.Z Ihsan wrote:
>
> Hi,
>
> i have question regarding the kdtree class in deal ii, though it has 
> become a deprecated class, but for my particular problem it might be 
> helpful. 
> suppose i have a class
>
> template
>
> class A
> {
> public: 
> A();
> 
> private:
> KDTree kdtree;
> };
>
>
>
> the question is how we can get a get function to the return the private 
> kdtree, e.g. A::get_kdtree() ? i still could not figure out what is the 
> type of the KDTree object. 
>
> BR, 
> Ihsan
>

-- 
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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/6f10ec6e-08ea-4fcd-aeb6-5db841d4b228%40googlegroups.com.

Re: [deal.II] Segfault when calling gather_evaluate() or read_dof_values_plain()

['Maxi Miller' via deal.II User Group]

For approximating the inverse mass matrix I followed step-48 with some 
modifications:
In the reinit-subroutine I added 
data.initialize_dof_vector(inv_mass_matrix);
FEEvaluation fe_eval(data);
const unsigned int   n_q_points = fe_eval.n_q_points;
for (unsigned int cell = 0; cell < data.n_cell_batches(); ++cell)
{
fe_eval.reinit(cell);
for (unsigned int q = 0; q < n_q_points; ++q)
fe_eval.submit_value(make_vectorized_array(1.), q);
fe_eval.integrate(true, false);
fe_eval.distribute_local_to_global(inv_mass_matrix);
}
inv_mass_matrix.compress(VectorOperation::add);
for (unsigned int k = 0; k < inv_mass_matrix.local_size(); ++k)
if (inv_mass_matrix.local_element(k) > 1e-15)
inv_mass_matrix.local_element(k) =
1. / inv_mass_matrix.local_element(k);
else
inv_mass_matrix.local_element(k) = 1;

and for applying the matrix I wrote
template 
void LaplaceOperator::local_apply_inverse_mass_matrix(
const MatrixFree &   data,
LinearAlgebra::distributed::Vector &  dst,
const LinearAlgebra::distributed::Vector ,
const std::pair & cell_range) 
const
{
dst.reinit(src);
dst = src;
dst.scale(inv_mass_matrix);
}
so that I could reuse the cell_loops from the original code. Still, the 
results indicate that something goes wrong. Should I handle the inverse 
matrix somehow different?
Thanks!
Am Dienstag, 14. Januar 2020 12:43:40 UTC+1 schrieb Martin Kronbichler:
>
> Dear Maxi,
>
> It looks like you are using a scalar finite element (FE_Q), but you 
> set FEEvaluation to operate on a vectorial field, i.e., 
>
> FEEvaluation phi(data);
>
> Notice the 4th template parameter "dim", which should be one. I agree it 
> is unfortunate that we do not provide a better error message, so I opened 
> an issue for it, https://github.com/dealii/dealii/issues/9312
>
> If you switch the parameter to 1, it should work. However, I want to point 
> out that you use a continuous element with CellwiseInverseMassMatrix - this 
> will not give you a valid matrix inverse. You need to either use a diagonal 
> mass matrix (see step-48) or use the cellwise inverse as a preconditioner 
> (when combined with suitable scaling) and solve the mass matrix iteratively.
>
> Best,
> Martin
> On 14.01.20 10:49, 'Maxi Miller' via deal.II User Group wrote:
>
> I could narrow it down to the function
>   
>  // same as above for has_partitioners_are_compatible == true
>   template <
> typename VectorType,
> typename std::enable_if::
> value,
> VectorType>::type * = nullptr>
>   inline void
>   check_vector_compatibility(
> const VectorType &vec,
> const internal::MatrixFreeFunctions::DoFInfo _info)
>   {
> (void)vec;
> (void)dof_info;
> Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
>ExcMessage(
>  "The parallel layout of the given vector is not "
>  "compatible with the parallel partitioning in MatrixFree. "
>  "Use MatrixFree::initialize_dof_vector to get a "
>  "compatible vector."));
>   }
>
> located in vector_access_internal.h. Apparently I have a segfault in the 
> function 
> 
>  template 
> inline bool
> Vector::partitioners_are_compatible(
>   const Utilities::MPI::Partitioner ) const
> {
>   return partitioner->is_compatible(part);
> }
>
> located in la_parallel_vector.templates.h, and thereby directly stopping 
> the program without triggering the assert()-function. The direct gdb output 
> for this function is 
> #3  0x700c0c5c in 
> dealii::LinearAlgebra::distributed::Vector dealii::MemorySpace::Host>::partitioners_are_compatible (this=0x0, 
> part=...) at 
> ~Downloads/git-files/dealii/include/deal.II/lac/la_parallel_vector.templates.h:2021
>
> I am not sure what that means. Could it point to a nullpointer somewhere 
> (after this=0x0)? Though, when putting a breakpoint to the first function 
> and printing the involved variables there, I get
> (gdb) print dof_info.vector_partitioner
> $2 = std::shared_ptr (use 
> count 3, weak count 0) = {get() = 0x7a60c0}
> (gdb) print vec.partitioner
> $3 = std::shared_ptr (use 
> count 3, weak count 0) = {get() = 0x7a60c0}
>
> i.e. no null ptr. Could the problem be there, or somewhere else?
> Thanks!
>
> Am Montag, 13. Januar 2020 21:41:54 UTC+1 schrieb Wolfgang Bangerth: 
>>
>> On 1/13/20 8:41 AM, 'Maxi Miller' via deal.II User Group wrote: 
>> > Therefore, I assume that I have some bugs in my code, but am not 
>> experienced 
>> > enough in writing matrix-free code for being able to debug it myself. 
>> What 
>> > kind of approach could I do now, or is there simply a glaring bug which 

Re: [deal.II] discontinous contour over elements

[Bruno Blais]

Dear David,
How are you calculating the vorticity?
As Wolfgang and Praveen have mentioned, if you are using the 
DataPostProcessor, then this will use your shape functions to calculate the 
vorticity. However, your P2-P1 elements are only C0 continuous. 
Consequently, your vorticity can possibly be inherently discontinuous at 
the element edges. However, I am surprised that you obtain such strong 
discontinuity. In our code based on deal.ii (
https://github.com/lethe-cfd/lethe) we implement stabilized formulations 
for the NS equations and the vorticity results for such cases are very 
smooth (even when represented using discontinuous shape functions.

Have you established the convergence of your code using manufactured 
solution? This is the first thing I would suggest. You can look at the 
applications_tests of lethe for some examples of easy manufactured 
solutions for the Incompressible Navier-Stokes equations. There are also 
common books that treat this issue (for instance : 
https://www.amazon.ca/Verification-Validation-Scientific-Computing-Oberkampf/dp/0521113601
)

Please feel free to reach out if you have more questions.
Best
Bruno




On Monday, 13 January 2020 20:29:53 UTC-5, David Eaton wrote:
>
> Thanks Wolfgang and Praveen for providing suggestions. I have tried to 
> debugging this code for a while. I have attached this simple code on this 
> email. I followed the instructions in tutorial closely. Hopefully, anyone 
> could give some suggestions.
>
> Best
> D. 
> --
> *From:* dea...@googlegroups.com   > on behalf of Wolfgang Bangerth  >
> *Sent:* Tuesday, January 14, 2020 6:24 AM
> *To:* dea...@googlegroups.com   >
> *Subject:* Re: [deal.II] discontinous contour over elements 
>  
> On 1/12/20 9:17 PM, David Eaton wrote:
> > My inflow condition is uniform. This formulation and mesh is tested in a 
> > simple C++ code without library. The  large mesh near the inflow does 
> not give 
> > this problem.
> > Yes. I am using C0 element. I did calculation using tecplot. However, 
> the 
> > results from a my C++ code does not give this issue either. Just now, I 
> check 
> > the formulation again. Although I use Q2Q1 Taylor-Hood element without 
> any 
> > stabilization, these issues are still happening.
>
> David -- we don't really know what formulation you are using, how you are 
> implementing it, what you are comparing against, and a number of other 
> factors.
>
> If you have a formulation that computes u,p, and you are plotting the 
> vorticity, you need to expect that the isocontours are discontinuous for 
> the 
> reasons Praveen already stated. If you are getting results that make no 
> sense 
> to you, then the first step would be to ensure that your program is 
> converging 
> as expected. To do this, choose a solution that you know and compute the 
> error 
> norm; then ensure that the program yields error norms that decrease as 
> expected with mesh refinement.
>
> Best
>   W.
>
> -- 
> 
> Wolfgang Bangerth  email: bang...@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 dea...@googlegroups.com .
> To view this discussion on the web visit 
> https://groups.google.com/d/msgid/dealii/e1b55e95-a139-d62e-e786-61393cf84dc7%40colostate.edu
> .
>

-- 
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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/5843231a-5b39-4f47-8b8e-ca3061959295%40googlegroups.com.

[deal.II] KDTree question

[A.Z Ihsan]

Hi,

i have question regarding the kdtree class in deal ii, though it has become 
a deprecated class, but for my particular problem it might be helpful. 
suppose i have a class

template

class A
{
public: 
A();

private:
KDTree kdtree;
};



the question is how we can get a get function to the return the private 
kdtree, e.g. A::get_kdtree() ? i still could not figure out what is the 
type of the KDTree object. 

BR, 
Ihsan

-- 
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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/b2f44c8e-5ff6-4a18-b382-e897e1f883bf%40googlegroups.com.

Re: [deal.II] Segfault when calling gather_evaluate() or read_dof_values_plain()

[Martin Kronbichler]

Dear Maxi,

It looks like you are using a scalar finite element (FE_Q), but you 
set FEEvaluation to operate on a vectorial field, i.e.,


        FEEvaluation phi(data);

Notice the 4th template parameter "dim", which should be one. I agree it 
is unfortunate that we do not provide a better error message, so I 
opened an issue for it, https://github.com/dealii/dealii/issues/9312


If you switch the parameter to 1, it should work. However, I want to 
point out that you use a continuous element with 
CellwiseInverseMassMatrix - this will not give you a valid matrix 
inverse. You need to either use a diagonal mass matrix (see step-48) or 
use the cellwise inverse as a preconditioner (when combined with 
suitable scaling) and solve the mass matrix iteratively.


Best,
Martin

On 14.01.20 10:49, 'Maxi Miller' via deal.II User Group wrote:

I could narrow it down to the function
|
// same as above for has_partitioners_are_compatible == true
template<
typenameVectorType,
typenamestd::enable_if::value,
VectorType>::type *=nullptr>
inlinevoid
  check_vector_compatibility(
constVectorType&                     vec,
constinternal::MatrixFreeFunctions::DoFInfo_info)
{
(void)vec;
(void)dof_info;
Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
ExcMessage(
"The parallel layout of the given vector is not "
"compatible with the parallel partitioning in MatrixFree. "
"Use MatrixFree::initialize_dof_vector to get a "
"compatible vector."));
}
|

located in vector_access_internal.h. Apparently I have a segfault in 
the function

|
template
inlinebool
Vector::partitioners_are_compatible(
constUtilities::MPI::Partitioner)const
{
returnpartitioner->is_compatible(part);
}
|

located in la_parallel_vector.templates.h, and thereby directly 
stopping the program without triggering the assert()-function. The 
direct gdb output for this function is

|
#3  0x700c0c5c in 
dealii::LinearAlgebra::distributed::Vectordealii::MemorySpace::Host>::partitioners_are_compatible (this=0x0, 
part=...) at 
~Downloads/git-files/dealii/include/deal.II/lac/la_parallel_vector.templates.h:2021

|

I am not sure what that means. Could it point to a nullpointer 
somewhere (after this=0x0)? Though, when putting a breakpoint to the 
first function and printing the involved variables there, I get

|
(gdb)printdof_info.vector_partitioner
$2 =std::shared_ptr(usecount 
3,weak count 0)={get()=0x7a60c0}

(gdb)printvec.partitioner
$3 =std::shared_ptr(usecount 
3,weak count 0)={get()=0x7a60c0}

|

i.e. no null ptr. Could the problem be there, or somewhere else?
Thanks!

Am Montag, 13. Januar 2020 21:41:54 UTC+1 schrieb Wolfgang Bangerth:

On 1/13/20 8:41 AM, 'Maxi Miller' via deal.II User Group wrote:
> Therefore, I assume that I have some bugs in my code, but am not
experienced
> enough in writing matrix-free code for being able to debug it
myself. What
> kind of approach could I do now, or is there simply a glaring
bug which I did
> not see yet?

I haven't even looked at the code yet, but for segfaults there is
a relatively
simple cure: Run your code in a debugger. A segmentation fault is
a problem
where some piece of code tries to access data at an address
(typically through
a pointer) that is not readable or writable. The general solution
to finding
out what the issue is is to run the code in a debugger -- the
debugger will
then stop at the place where the problem happens, and you can
inspect all of
the values of the pointers and variables that are live at that
location. Once
you see what variable is at fault, the problem is either already
obvious, or
you can start to think about how a pointer got the value it has
and why.

Best
  W.

-- 



Wolfgang Bangerth          email: bang...@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 
.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/d55e53b2-a0c8-4cc0-ab84-7521f8b23d81%40googlegroups.com 
.


--
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 

Re: [deal.II] Segfault when calling gather_evaluate() or read_dof_values_plain()

['Maxi Miller' via deal.II User Group]

I could narrow it down to the function
 
 // same as above for has_partitioners_are_compatible == true
  template <
typename VectorType,
typename std::enable_if::
value,
VectorType>::type * = nullptr>
  inline void
  check_vector_compatibility(
const VectorType &vec,
const internal::MatrixFreeFunctions::DoFInfo _info)
  {
(void)vec;
(void)dof_info;
Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
   ExcMessage(
 "The parallel layout of the given vector is not "
 "compatible with the parallel partitioning in MatrixFree. "
 "Use MatrixFree::initialize_dof_vector to get a "
 "compatible vector."));
  }

located in vector_access_internal.h. Apparently I have a segfault in the 
function 
   
 template 
inline bool
Vector::partitioners_are_compatible(
  const Utilities::MPI::Partitioner ) const
{
  return partitioner->is_compatible(part);
}

located in la_parallel_vector.templates.h, and thereby directly stopping 
the program without triggering the assert()-function. The direct gdb output 
for this function is 
#3  0x700c0c5c in 
dealii::LinearAlgebra::distributed::Vector::partitioners_are_compatible (this=0x0, 
part=...) at 
~Downloads/git-files/dealii/include/deal.II/lac/la_parallel_vector.templates.h:2021

I am not sure what that means. Could it point to a nullpointer somewhere 
(after this=0x0)? Though, when putting a breakpoint to the first function 
and printing the involved variables there, I get
(gdb) print dof_info.vector_partitioner
$2 = std::shared_ptr (use count 3
, weak count 0) = {get() = 0x7a60c0}
(gdb) print vec.partitioner
$3 = std::shared_ptr (use count 3
, weak count 0) = {get() = 0x7a60c0}

i.e. no null ptr. Could the problem be there, or somewhere else?
Thanks!

Am Montag, 13. Januar 2020 21:41:54 UTC+1 schrieb Wolfgang Bangerth:
>
> On 1/13/20 8:41 AM, 'Maxi Miller' via deal.II User Group wrote: 
> > Therefore, I assume that I have some bugs in my code, but am not 
> experienced 
> > enough in writing matrix-free code for being able to debug it myself. 
> What 
> > kind of approach could I do now, or is there simply a glaring bug which 
> I did 
> > not see yet? 
>
> I haven't even looked at the code yet, but for segfaults there is a 
> relatively 
> simple cure: Run your code in a debugger. A segmentation fault is a 
> problem 
> where some piece of code tries to access data at an address (typically 
> through 
> a pointer) that is not readable or writable. The general solution to 
> finding 
> out what the issue is is to run the code in a debugger -- the debugger 
> will 
> then stop at the place where the problem happens, and you can inspect all 
> of 
> the values of the pointers and variables that are live at that location. 
> Once 
> you see what variable is at fault, the problem is either already obvious, 
> or 
> you can start to think about how a pointer got the value it has and why. 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/d55e53b2-a0c8-4cc0-ab84-7521f8b23d81%40googlegroups.com.

Re: [deal.II] Using solution on from one FE problem as a boundary condition for another

[Ernesto Ismail]

Dear Prof. Bangerth,

Thank you for your reply. 

I would ideally have liked to impose the boundary conditions strongly but 
it looks like doing so will not be straight-forward. I will do some 
investigations into applying the boundary conditions weakly and attempt 
that first.

One approach I had tried was to use Lobatto quadrature and use a quadrature 
point history data structure to transfer the solution at the boundary. The 
problem with this approach was that the differing element types lead to 
some strange behaviours along the boundary.

Thanks again,
Ernesto

On Thursday, 9 January 2020 23:43:27 UTC+2, Wolfgang Bangerth wrote:
>
>
> Ernesto, 
>
> > I am currently working with a staggered coupled system where I am 
> solving 
> > several finite element problems in sequence. Up until now my boundary 
> > conditions have been fairly straightforward and I've been able to get 
> > meaningful results from my work. 
> > 
> > The finite element schemes I use for each step are fairly different, 
> with some 
> > being discontinuous and of varying orders. 
> > 
> > I am now in a position where I need to use the solution of one of my 
> steps as 
> > the boundary condition for the next (I actually need to combine two of 
> them 
> > with a simple arithmetic function, but baby steps for now). I presume I 
> need 
> > to use some sort of projection to achieve this, but I'm not exactly sure 
> how 
> > to go about it. I thought to try and use the project_boundary_values() 
> > function to do what I need, but I am unsure how I would structure 
> > the boundary_functions argument. 
>
> Do you need to enforce these boundary conditions strongly? If you can 
> enforce 
> them weakly, then the solution of one problem only enters the weak form of 
> another problem (in that case through a boundary integral) and that is 
> really 
> not very different from the way we do this in other problems where we have 
> multiple DoFHandlers -- e.g., in step-31. (Though my usual advice still 
> holds 
> true that I think everything becomes simpler if you pack everything into 
> one 
> DoFHandler). 
>
> It gets more complicated if you want to impose strong boundary conditions. 
> You've already found the FEFieldFunction function, but it's complicated 
> and 
> slow. A simpler way would be if you looked at the implementation of the 
> interpolate_boundary_values() function and identified where it is asking 
> the 
> Function object for its values. That's the place where you'd need to 
> somehow 
> combine your existing solutions into something else and turn that into a 
> constraint. 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/832bb83a-7e04-4243-9524-64db1bfc9345%40googlegroups.com.