Re: [deal.II] The shape function on physical cell in hexahedral is not parallel with the edge (FeNedelecSZ and FENedelec)

2019-04-15 Thread Phạm Ngọc Kiên
Hi,
I have written this code for testing purpose.
In case of a rotated cube, the shape functions seemed to be good as they
were parallel to the edge, although there was some small round-off errors
at the value 0.
However, when I tested with the grid which was loaded from GMSH, the shape
functions failed.
In the attachment, I showed the information for the first active cell only.
It is easy to get information of other cells with the codes.

I think that there exists something wrong with the cell that are not a
cube, but a hexahedron when initializing the shape functions (fe_values).

I think it would be easier for you to help me with the attachment below.
Thank you very much.
Best regards.
Pham Ngoc Kien

Vào Th 2, 15 thg 4, 2019 vào lúc 21:48 Wolfgang Bangerth <
bange...@colostate.edu> đã viết:

> On 4/11/19 9:17 PM, Phạm Ngọc Kiên wrote:
> > Testing for an edge whose global vertices located from (0,0,0) to
> (0,0,1) in
> > real coordinates.
> > With a cube I get the shape function vectors at the dof related to the
> edge,
> > for examples, (0,0,0), (0,0,-0.25), (0,0,-0.5), (0,0,-1), which are
> parallel
> > to the edge.
> > However, if I create a mesh with GMSH and check the shape function
> vectors,
> > some of them are: (-0.0868383, -0.428, 1), (-0.129933, -0.402772, 1),
> > (-0.052103 -0.2568 0.6),..,
> > which are not parallel to the edge.
> >
> > In my limited opinion, I will get the wrong solution as the shape
> functions
> > created by my code are not parallel to the edge in my model.
> > Could you please tell me why this happen and how to fix it?
>
> Good question. This does sound wrong. To figure out what is happening
> exactly,
> we typically construct small test cases that try to illustrate one
> particular
> issue. In your case, I would try to create a mesh with exactly one cell --
> either the one you have above, or maybe even simpler a cube that is just
> rotated. Then create a DoFHandler on it and output the values of the shape
> functions. This way, you really have only one thing that could go wrong,
> and
> it's easy to demonstrate what you see and how it differs from expectations.
>
> Do you think you could come up with a small program that does this?
>
> Best
>   WB
>
> --
> 
> 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.
>

-- 
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.
/* -
 *
 * Copyright (C) 2013 - 2018 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE.md at
 * the top level directory of deal.II.
 *
 * -

 *
 * Author: Pham Ngoc Kien, Seoul National University, 2019
 */


#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 


#include 
#include 
#include 
#include 

#include 
#include 
#include 
#include 


#include 
#include 
#include 
#include 
#include 

#include 
#include 
#include 
#include 


#include 
#include 

#include 


#include 
#include 
#include 

#include 
#include 
#include 

using namespace dealii;

template
class EM {
public:
EM();

~EM();

void run();


private:
void grid_generation();

void setup_system();

void assemble_system();


void print_mesh_info(const Triangulation ,
 const std::string );


Triangulation triangulation;
MappingQ mapping;
const unsigned int order = 0; //order p of element
DoFHandler dof_handler;
FESystem fe;

Pointget_vertex_pos(const typename DoFHandler::active_cell_iterator ,
 const FESystem 

Re: [deal.II] step-22 compute_no_normal_flux_constraints

2019-04-15 Thread Timo Heister
For a weak solution, you can split the boundary integral (coming from
integration by parts) into a normal and a tangential part. With
no-normal flux contraints you remove the normal part only, so the
tangential part needs to be written as a boundary integral on the RHS
of your PDE unless it is zero in your manufactured solution.
You can find more info in Ern/Guermond in the "variations on boundary
conditions" for Stokes.

On Mon, Apr 15, 2019 at 10:01 AM  wrote:
>
> Hi Wolfgang
>
> Thanks for your reply.
>
> That was actually what I had done previously. I've tried all sorts by looking 
> at the values on different boundaries, indicators, tried multiple problems, 
> etc before I had asked the question, which is why I came to thinking it had 
> something with the implementation.
>
> Thanks for your suggestions - I had another thought, that perhaps it might be 
> a mathematical issue. On a no normal flux condition in the Stokes case, am I 
> correct in assuming that for the reamining components (1- nxn)(n.[pI-2e(u)]) 
> for the boundary ocndition, if homogenous (=0) then it disappears within the 
> weak form, and if non zero, it needs to be imposed weakly?
>
> or does something else need to take place for this?
>
> On Monday, April 15, 2019 at 3:33:36 PM UTC+1, Wolfgang Bangerth wrote:
>>
>>
>> > The mesh has over 10 cells in it, it is super refined. And oddly, 
>> > when
>> > the refinement level is less, it doesn't blow up. It's only after a
>> > certain point.
>> > it is an even global refinement, starting from a hyper divided 
>> > rectangle.
>> > no fancy refinement.
>> >
>> > I only realised this as I wasn't getting the right convergence rates 
>> > for v
>> > at the lower refinement levels so I just kept going up to see if it
>> > eventually did and it blew up.
>> > Until this point, at the refinement levels that work, the pressure
>> > converges correctly at two consecutive degrees, and so does v at the
>> > higher degree. but at the lower degree, the convergence rates for v 
>> > were
>> > decreasing with refinement.
>> > That;s why I continued to refine then came across this situation.
>> >
>> > Seeing that it is only at the bottom, and noticing that if I used
>> > Dirichlet condition on that boundary, then the solution doesn't blow up
>> > and converges correctly, I can only assume it's to do with the normal 
>> > flux
>> > condition there, which i think I'm imposing correctly.
>>
>> Jane -- I don't think anyone other than you is in a position to figure this
>> out. I've learned to debug problems with a mind set of assuming that
>> everything I believe is correct may in fact not. So when you say
>>"I can only assume it's to do with the normal flux condition there, which
>> i think I'm imposing correctly"
>> then my approach is to assume that they are not imposed correctly, and to
>> write code that helps me to *verify* that it is.
>>
>> As an example, even though you think that you've set the boundary indicators
>> correctly at all cells on the bottom, this may not be the case for whatever
>> reason. Write a loop over all cells and all boundary faces and *output*
>> location and boundary indicator, to make sure it is. Go through this sort of
>> procedure for everything you *believe* should be true to make sure that it
>> *is* true.
>>
>> I'm sorry I have no other suggestions for what to do. We've all been in your
>> situation where stuff doesn't work for reasons that seem completely
>> mysterious. The successful among us are the ones who have built the mental 
>> and
>> computer skills to figure out what the cause is. Among the mental skills is
>> the ability to assume that what one believes to be true (because one has
>> written the code) may not be so. The computer skills is then to use a 
>> debugger
>> or print statements in the right places to verify that something is indeed 
>> the
>> case or not, and to narrow down the possible root causes.
>>
>> Best
>>   W.
>>
>>
>> --
>> 
>> Wolfgang Bangerth  email: bang...@colostate.edu
>> www: 
>> https://urldefense.proofpoint.com/v2/url?u=http-3A__www.math.colostate.edu_-7Ebangerth_=DwIFaQ=Ngd-ta5yRYsqeUsEDgxhcqsYYY1Xs5ogLxWPA_2Wlc4=4k7iKXbjGC8LfYxVJJXiaYVu6FRWmEjX38S7JmlS9Vw=QVz-0qa-BE17tx8qPdxGM1rRo9ff7D68_fU_ptB7amo=54A9voRgVztKomopsSD34VYdFjTc_5ZY7Lj3-cyi2SE=
>>
> --
> The deal.II project is located at 
> https://urldefense.proofpoint.com/v2/url?u=http-3A__www.dealii.org_=DwIFaQ=Ngd-ta5yRYsqeUsEDgxhcqsYYY1Xs5ogLxWPA_2Wlc4=4k7iKXbjGC8LfYxVJJXiaYVu6FRWmEjX38S7JmlS9Vw=QVz-0qa-BE17tx8qPdxGM1rRo9ff7D68_fU_ptB7amo=2YX98RoE1A9bLcn1QbLuBECzNPFA4vb6_oBwLDbzKR0=
> For mailing list/forum options, see 
> 

Re: [deal.II] step-22 compute_no_normal_flux_constraints

2019-04-15 Thread jane . lee
Hi Wolfgang 

Thanks for your reply. 

That was actually what I had done previously. I've tried all sorts by 
looking at the values on different boundaries, indicators, tried multiple 
problems, etc before I had asked the question, which is why I came to 
thinking it had something with the implementation.

Thanks for your suggestions - I had another thought, that perhaps it might 
be a mathematical issue. On a no normal flux condition in the Stokes case, 
am I correct in assuming that for the reamining components (1- 
nxn)(n.[pI-2e(u)]) for the boundary ocndition, if homogenous (=0) then it 
disappears within the weak form, and if non zero, it needs to be imposed 
weakly?

or does something else need to take place for this? 

On Monday, April 15, 2019 at 3:33:36 PM UTC+1, Wolfgang Bangerth wrote:
>
>
> > The mesh has over 10 cells in it, it is super refined. And 
> oddly, when 
> > the refinement level is less, it doesn't blow up. It's only after a 
> > certain point. 
> > it is an even global refinement, starting from a hyper divided 
> rectangle. 
> > no fancy refinement. 
> > 
> > I only realised this as I wasn't getting the right convergence rates 
> for v 
> > at the lower refinement levels so I just kept going up to see if it 
> > eventually did and it blew up. 
> > Until this point, at the refinement levels that work, the pressure 
> > converges correctly at two consecutive degrees, and so does v at the 
> > higher degree. but at the lower degree, the convergence rates for v 
> were 
> > decreasing with refinement. 
> > That;s why I continued to refine then came across this situation. 
> > 
> > Seeing that it is only at the bottom, and noticing that if I used 
> > Dirichlet condition on that boundary, then the solution doesn't blow 
> up 
> > and converges correctly, I can only assume it's to do with the 
> normal flux 
> > condition there, which i think I'm imposing correctly. 
>
> Jane -- I don't think anyone other than you is in a position to figure 
> this 
> out. I've learned to debug problems with a mind set of assuming that 
> everything I believe is correct may in fact not. So when you say 
>"I can only assume it's to do with the normal flux condition there, 
> which 
> i think I'm imposing correctly" 
> then my approach is to assume that they are not imposed correctly, and to 
> write code that helps me to *verify* that it is. 
>
> As an example, even though you think that you've set the boundary 
> indicators 
> correctly at all cells on the bottom, this may not be the case for 
> whatever 
> reason. Write a loop over all cells and all boundary faces and *output* 
> location and boundary indicator, to make sure it is. Go through this sort 
> of 
> procedure for everything you *believe* should be true to make sure that it 
> *is* true. 
>
> I'm sorry I have no other suggestions for what to do. We've all been in 
> your 
> situation where stuff doesn't work for reasons that seem completely 
> mysterious. The successful among us are the ones who have built the mental 
> and 
> computer skills to figure out what the cause is. Among the mental skills 
> is 
> the ability to assume that what one believes to be true (because one has 
> written the code) may not be so. The computer skills is then to use a 
> debugger 
> or print statements in the right places to verify that something is indeed 
> the 
> case or not, and to narrow down the possible root causes. 
>
> 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.
For more options, visit https://groups.google.com/d/optout.


Re: [deal.II] Accuracy of Dirichlet condition for p in step-20

2019-04-15 Thread jane . lee
Hi Wolfgang,

I can get a table if it would be useful. 

I see what you mean in terms of convergence. I guess I was looking for the 
accuracy pointwise on a boundary where the Dirichlet condition for the 
pressure is imposed weakly. In my case, the value of the output on the 
boundary was important, which is why I was looking at the value of the 
pressure on the boundary. 

I will try using Q2-Q2 elements, though I realise that this doesn't change 
the application of the BC itself, so the issues may be the same, but 
definitely worth a try. 

For a neumann boundary condition though (my problem has that on another 
boundary), am I correct in thinking that project_boundary_values doesn't 
allow for component masking? the condition is only on u = f(grad p) and not 
on p itself so I would need something similar to component masking the 
velocities. 

Many thanks
 

On Monday, April 15, 2019 at 3:23:28 PM UTC+1, Wolfgang Bangerth wrote:
>
>
> Jane, 
>
> > I continued to find out why I wasn't getting the correct applied 
> Dirichlet 
> > values on the boundary for a code very similar to step-20, where the 
> Dirichlet 
> > condition is applied weakly using 
> > 
> > for (unsigned int face_no=0; 
> > face_no::faces_per_cell; 
> > ++face_no) 
> > if (cell->at_boundary(face_no)) 
> > { 
> > fe_face_values.reinit 
> > <
> https://www.dealii.org/8.4.1/doxygen/deal.II/classFEFaceValues.html#af6e079ca7429d54433343d50bd334c3c>
>  
>
> > (cell, face_no); 
> > pressure_boundary_values 
> > .value_list (fe_face_values.get_quadrature_points 
> > <
> https://www.dealii.org/8.4.1/doxygen/deal.II/classFEValuesBase.html#a5f8732ebe2d3c6746f6de26a79cb1e45>(),
>  
>
> > boundary_values); 
> > for (unsigned int q=0; q > for (unsigned int i=0; i > local_rhs(i) += -(fe_face_values[velocities].value (i, q) * 
> > fe_face_values.normal_vector 
> > <
> https://www.dealii.org/8.4.1/doxygen/deal.II/classFEValuesBase.html#a130eea0fa89263d93b20521addc830c7>(q)
>  
>
> > * 
> > boundary_values[q] * 
> > fe_face_values.JxW 
> > <
> https://www.dealii.org/8.4.1/doxygen/deal.II/classFEValuesBase.html#ad097580a2f71878695096cc73b271b9d>(q));
>  
>
> > } 
> > 
> > 
> > I then looked at step-20 - I used the exact code but solved directly 
> instead, 
> > giving me the same results as in the tutorial (for the errors etc). 
> > 
> > Even in step-20, the boundary values aren't correct for most of the 
> tests, or 
> > rather, have a lot of error in itself. 
> > For example, at the point (1,1) which is on the boundary, p = -1.1 with 
> the 
> > given test problem. 
> > 
> > These are the values I obtained having run the code: 
> > at lowest degree (0), refinement level (RL) 3, p=-0.941992, which is 
> rather 
> > far off the -1.1 value for the Dirichlet condition applied. even at RL6, 
> p=1.07976 
> > I tried the next degree up (1), at RL3, p= -1.09984. eventually at RL6 
> for 
> > this degree, we get -1.1. 
> > Similarly, at degree 2, at RL3. p is still not accurate at p = -1.10004 
>
> Individual values are usually not particularly meaningful. What matters is 
> the 
> convergence towards the correct value. Can you produce a table for each 
> order 
> in which you show the value as a function of refinement level? 
>
> I will add that it's not clear to me that the finite element method 
> applied to 
> this problem guarantees that the error *at a single point* converges to 
> zero. 
> We know that the error converges to zero when measured in integral 
> quantities 
> such as the L2 norm, but that mean not imply pointwise convergence. 
>
>
> > How else can we imposed the Dirichlet condition (without having to use 
> very 
> > high refinement levels and degrees) on the boundary for this problem? 
>
> I don't know of other ways. The question I would ask first is whether what 
> you 
> have *converges*. That's all you can really hope for. The second question 
> is 
> whether the error is *small*. If the method converges, then the error 
> *will* 
> be small if only the mesh is fine enough; but "fine enough" may be smaller 
> than you can computationally afford, and in that case one can start to ask 
> about alternatives (such as using the Q2-Q1 element used in step-43 
> instead of 
> the RT element you have here), but the first step is to unambiguously 
> establish convergence. 
>
> Best 
>   WB 
>
> -- 
>  
> 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.
For more options, visit 

Re: [deal.II] Mass matrix for a distributed vector problem

2019-04-15 Thread Robert Spartus
Dear Wolfgang,

> It is hard to imagine situations in which the mass matrix would be
singular.
> It is a positive definite form that gives rise to the mass matrix and so
it
> really shouldn't be singular at all. Can you show the code again with
which
> you build it?

It seems that my mesh is neither singular nor degenerate. I wonder if the
problem is that I am solving a vector valued problem. To build this mass
matrix, I adapted the function Diffusion::assemble_system from step-52.

You will find the code attached, as well as the output of one run. If I
isolate the mass matrix built and calculate its determinant in Python, the
result is indeed zero.

I am really at a loss here, so any help is appreciated.

Bests,
Bob

On Mon, 15 Apr 2019 at 16:26, Wolfgang Bangerth 
wrote:

> On 4/14/19 11:59 PM, Robert Spartus wrote:
> >
> > Thanks for the insightful discussion on the integrating issue. Wolfgang,
> I
> > guess your last argument is the same as you gave in one of your
> fantastic
> > lectures?
>
> Yes. (Also, thanks for the compliment :-) )
>
>
> > Incidentally, do you have any ideas on how to solve the singularity of
> the
> > mass matrix in this vector-valued problem?
>
> It is hard to imagine situations in which the mass matrix would be
> singular.
> It is a positive definite form that gives rise to the mass matrix and so
> it
> really shouldn't be singular at all. Can you show the code again with
> which
> you build it?
>
> The only situation where the mass matrix may become singular is if you
> have a
> degenerate mesh with elements of zero or negative volume.
>
> Best
>   WB
>
>
> --
> 
> 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.
>

-- 
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.
 2.8e-02  0.0e+00  0.0e+00  0.0e+00  0.0e+00  0.0e+00  
0.0e+00  0.0e+00  0.0e+00  0.0e+00  0.0e+00  0.0e+00  
0.0e+00  0.0e+00  0.0e+00  0.0e+00  0.0e+00  0.0e+00 
 0.0e+00  2.8e-02  0.0e+00  0.0e+00  0.0e+00  0.0e+00  
0.0e+00  0.0e+00  0.0e+00  0.0e+00  0.0e+00  0.0e+00  
0.0e+00  0.0e+00  0.0e+00  0.0e+00  0.0e+00  0.0e+00 
 0.0e+00  0.0e+00  5.6e-02  0.0e+00  0.0e+00  6.9e-03  
2.8e-02  2.8e-02  0.0e+00  0.0e+00  0.0e+00  6.9e-03  
0.0e+00  0.0e+00  0.0e+00  0.0e+00  0.0e+00  0.0e+00 
 0.0e+00  0.0e+00  0.0e+00  5.6e-02  0.0e+00  0.0e+00  
0.0e+00  0.0e+00  0.0e+00  0.0e+00  0.0e+00  0.0e+00  
0.0e+00  0.0e+00  0.0e+00  0.0e+00  0.0e+00  0.0e+00 
 0.0e+00  0.0e+00  0.0e+00  0.0e+00  5.6e-02  0.0e+00  
0.0e+00  0.0e+00  0.0e+00  0.0e+00  0.0e+00  0.0e+00  
0.0e+00  0.0e+00  0.0e+00  0.0e+00  0.0e+00  0.0e+00 
 0.0e+00  0.0e+00  6.9e-03  0.0e+00  0.0e+00  5.6e-02  
2.8e-02  2.8e-02  0.0e+00  0.0e+00  0.0e+00  0.0e+00  
0.0e+00  0.0e+00  6.9e-03  0.0e+00  0.0e+00  0.0e+00 
 0.0e+00  0.0e+00  2.8e-02  0.0e+00  0.0e+00  2.8e-02  
1.1e-01  1.1e-01  0.0e+00  0.0e+00  0.0e+00  2.8e-02  
0.0e+00  0.0e+00  2.8e-02  0.0e+00  0.0e+00  0.0e+00 
 0.0e+00  0.0e+00  2.8e-02  0.0e+00  0.0e+00  2.8e-02  
1.1e-01  1.1e-01  0.0e+00  0.0e+00  0.0e+00  2.8e-02  
0.0e+00  0.0e+00  2.8e-02  0.0e+00  0.0e+00  0.0e+00 
 0.0e+00  0.0e+00  0.0e+00  0.0e+00  0.0e+00  0.0e+00  
0.0e+00  0.0e+00  2.8e-02  0.0e+00  0.0e+00  0.0e+00  
0.0e+00  0.0e+00  0.0e+00  0.0e+00  0.0e+00  0.0e+00 
 0.0e+00  0.0e+00  0.0e+00  0.0e+00  0.0e+00  0.0e+00  
0.0e+00  0.0e+00  0.0e+00  2.8e-02  

Re: [deal.II] step-22 compute_no_normal_flux_constraints

2019-04-15 Thread Wolfgang Bangerth


> The mesh has over 10 cells in it, it is super refined. And oddly, when
> the refinement level is less, it doesn't blow up. It's only after a
> certain point.
> it is an even global refinement, starting from a hyper divided rectangle.
> no fancy refinement.
> 
> I only realised this as I wasn't getting the right convergence rates for v
> at the lower refinement levels so I just kept going up to see if it
> eventually did and it blew up.
> Until this point, at the refinement levels that work, the pressure
> converges correctly at two consecutive degrees, and so does v at the
> higher degree. but at the lower degree, the convergence rates for v were
> decreasing with refinement.
> That;s why I continued to refine then came across this situation.
> 
> Seeing that it is only at the bottom, and noticing that if I used
> Dirichlet condition on that boundary, then the solution doesn't blow up
> and converges correctly, I can only assume it's to do with the normal flux
> condition there, which i think I'm imposing correctly.

Jane -- I don't think anyone other than you is in a position to figure this 
out. I've learned to debug problems with a mind set of assuming that 
everything I believe is correct may in fact not. So when you say
   "I can only assume it's to do with the normal flux condition there, which
i think I'm imposing correctly"
then my approach is to assume that they are not imposed correctly, and to 
write code that helps me to *verify* that it is.

As an example, even though you think that you've set the boundary indicators 
correctly at all cells on the bottom, this may not be the case for whatever 
reason. Write a loop over all cells and all boundary faces and *output* 
location and boundary indicator, to make sure it is. Go through this sort of 
procedure for everything you *believe* should be true to make sure that it 
*is* true.

I'm sorry I have no other suggestions for what to do. We've all been in your 
situation where stuff doesn't work for reasons that seem completely 
mysterious. The successful among us are the ones who have built the mental and 
computer skills to figure out what the cause is. Among the mental skills is 
the ability to assume that what one believes to be true (because one has 
written the code) may not be so. The computer skills is then to use a debugger 
or print statements in the right places to verify that something is indeed the 
case or not, and to narrow down the possible root causes.

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] Mass matrix for a distributed vector problem

2019-04-15 Thread Wolfgang Bangerth
On 4/14/19 11:59 PM, Robert Spartus wrote:
> 
> Thanks for the insightful discussion on the integrating issue. Wolfgang, I 
> guess your last argument is the same as you gave in one of your fantastic 
> lectures?

Yes. (Also, thanks for the compliment :-) )


> Incidentally, do you have any ideas on how to solve the singularity of the 
> mass matrix in this vector-valued problem?

It is hard to imagine situations in which the mass matrix would be singular. 
It is a positive definite form that gives rise to the mass matrix and so it 
really shouldn't be singular at all. Can you show the code again with which 
you build it?

The only situation where the mass matrix may become singular is if you have a 
degenerate mesh with elements of zero or negative volume.

Best
  WB


-- 

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] Accuracy of Dirichlet condition for p in step-20

2019-04-15 Thread Wolfgang Bangerth


Jane,

> I continued to find out why I wasn't getting the correct applied Dirichlet 
> values on the boundary for a code very similar to step-20, where the 
> Dirichlet 
> condition is applied weakly using
> 
> 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);
> pressure_boundary_values
> .value_list (fe_face_values.get_quadrature_points 
> (),
> boundary_values);
> for (unsigned int q=0; q for (unsigned int i=0; i local_rhs(i) += -(fe_face_values[velocities].value (i, q) *
> fe_face_values.normal_vector 
> (q)
>  
> *
> boundary_values[q] *
> fe_face_values.JxW 
> (q));
> }
> 
> 
> I then looked at step-20 - I used the exact code but solved directly instead, 
> giving me the same results as in the tutorial (for the errors etc).
> 
> Even in step-20, the boundary values aren't correct for most of the tests, or 
> rather, have a lot of error in itself.
> For example, at the point (1,1) which is on the boundary, p = -1.1 with the 
> given test problem.
> 
> These are the values I obtained having run the code:
> at lowest degree (0), refinement level (RL) 3, p=-0.941992, which is rather 
> far off the -1.1 value for the Dirichlet condition applied. even at RL6, 
> p=1.07976
> I tried the next degree up (1), at RL3, p= -1.09984. eventually at RL6 for 
> this degree, we get -1.1.
> Similarly, at degree 2, at RL3. p is still not accurate at p = -1.10004

Individual values are usually not particularly meaningful. What matters is the 
convergence towards the correct value. Can you produce a table for each order 
in which you show the value as a function of refinement level?

I will add that it's not clear to me that the finite element method applied to 
this problem guarantees that the error *at a single point* converges to zero. 
We know that the error converges to zero when measured in integral quantities 
such as the L2 norm, but that mean not imply pointwise convergence.


> How else can we imposed the Dirichlet condition (without having to use very 
> high refinement levels and degrees) on the boundary for this problem?

I don't know of other ways. The question I would ask first is whether what you 
have *converges*. That's all you can really hope for. The second question is 
whether the error is *small*. If the method converges, then the error *will* 
be small if only the mesh is fine enough; but "fine enough" may be smaller 
than you can computationally afford, and in that case one can start to ask 
about alternatives (such as using the Q2-Q1 element used in step-43 instead of 
the RT element you have here), but the first step is to unambiguously 
establish convergence.

Best
  WB

-- 

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] The shape function on physical cell in hexahedral is not parallel with the edge (FeNedelecSZ and FENedelec)

2019-04-15 Thread Wolfgang Bangerth
On 4/11/19 9:17 PM, Phạm Ngọc Kiên wrote:
> Testing for an edge whose global vertices located from (0,0,0) to (0,0,1) in 
> real coordinates.
> With a cube I get the shape function vectors at the dof related to the edge, 
> for examples, (0,0,0), (0,0,-0.25), (0,0,-0.5), (0,0,-1), which are parallel 
> to the edge.
> However, if I create a mesh with GMSH and check the shape function vectors, 
> some of them are: (-0.0868383, -0.428, 1), (-0.129933, -0.402772, 1),  
> (-0.052103 -0.2568 0.6),..,
> which are not parallel to the edge.
> 
> In my limited opinion, I will get the wrong solution as the shape functions 
> created by my code are not parallel to the edge in my model.
> Could you please tell me why this happen and how to fix it?

Good question. This does sound wrong. To figure out what is happening exactly, 
we typically construct small test cases that try to illustrate one particular 
issue. In your case, I would try to create a mesh with exactly one cell -- 
either the one you have above, or maybe even simpler a cube that is just 
rotated. Then create a DoFHandler on it and output the values of the shape 
functions. This way, you really have only one thing that could go wrong, and 
it's easy to demonstrate what you see and how it differs from expectations.

Do you think you could come up with a small program that does this?

Best
  WB

-- 

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] Computation of Eigenvalues by Power Method with deal.ii

2019-04-15 Thread Wolfgang Bangerth
On 4/14/19 4:29 AM, illi wrote:
> I have the following code snippet for computing Eigenvalues using Power 
> Method:
> |
> Vector x;
> x = solution;
> double v = 0.0;
> PrimitiveVectorMemory> mem;
> const EigenPower>::AdditionalData data(0.);
> EigenPower<> ep(solver_control, mem, data);
> ep.solve(v, system_matrix, x);
> |
> (where system_matrix is a SparseMatrix)
> 
> The code arises an*abort trap*. What is wrong?

It's hard to tell without more information :-) Can you show the complete error 
message you see on the screen? I suspect you've already verified that the 
error really happens in the code snippet above, right? And that you are 
running in debug mode?

The right approach to debugging this will be to run the program in a debugger 
and see where the trap happens exactly.

Best
  WB


-- 

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.