[deal.II] make_sparsity_pattern, make_flux_sparsity_pattern problem when elements are connected with vertex

2019-03-19 Thread Yidong ZHAO
Hello all,

While using a new kind of shape function, degrees of freedom on each cell 
will couple to the degrees of freedom on other cells connected to the 
current one by a common face or common vertex.

>From the description of make_flux_sparsity_pattern(), it can solve the 
problem when two elements are connected by a common face, as depicted in 
elements_connected_by_common_line.png.

But if it's the case as elements_connected_by_common_vertex.png, is there 
any way to make_vertex_sparsity_pattern() ? 


Many thanks,
Yidong

-- 
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] How to build SparsityPattern when quadrature point has non-zero shape value to adjacent element's nodes

2019-03-19 Thread Yidong ZHAO
Hello,

While implementing a new kind shape function which has non-zero shape value 
between a quadrature point and an adjacent element's node, I found I 
couldn't find a way to make a SparsityPattern.

Could somebody give some suggestions how to build SparsityPattern?

The attach file is a picture shows that a quadrature point in an element 
will have non-zero shape value to an adjacent element's node.

Many thanks,
Yidong


-- 
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] Why the FENedelec and FENedelecSZ 's shape functions change versus edge length?

2019-03-19 Thread Phạm Ngọc Kiên
Hi colleagues,
When  testing my codes with FENedelec and FENedelecSZ elements, I saw that 
the shape functions change versus the mesh size.
I think that for each edge of a cell, the shape functions for the degree of 
freedom related to that edge are scaled with the inverse edge length.
For example, assuming I have a cell that is a cube, whose edge length is 1, 
the shape function at the quadrature point located on the edge is {0, 0, 1}.
For the smaller cell whose edge length is 0.5, the shape function values is 
{0, 0, 2}. 

In my limited opinion, the shape function is defined in unit cell over 
[0,1]^3, and the shape function should not change when the real cell change 
its size.
This problem leads to the change of my solution when using the finer mesh 
compared to that of coarser one. 
Because both my cell matrix and cell rhs change versus the edge length. 

For FeNedelecSZ, I tested with a cell limited in the region [0,1]^3, which 
means that the edge length is 1 everywhere.
When printing out the integral over the cell for one degree of freedom by 
   
\int_Omega  \varphi_i(x)  Js(x)  dx 
For simplicity, I set Omega =1 and Js(x) = {1,1,1} for testing the integral 
over the cell at degree of freedom i_th when source magnitude is 1.
The result for this is 0.25 rather than 1.
Is this means that for a source of magnitude 1., I actually can model it as 
a source edge at degree of freedom i_th with the magnitude 0.25?

I think there should be some idea that we change the shape function values 
like this.
Could you please give me some advice?

Thank you very much for your great support.
Best regards,

-- 
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] advection diffusion with periodic boundaries

2019-03-19 Thread Gary Uppal
Hi,

I'm trying to solve a time dependent advection-diffusion equation with 
periodic boundary conditions. Just a simple du/dt = D \nabla^2 u - v \dot 
\grad u for now. I use a large diffusion constant, so stability shouldn't 
be an issue. The solution behaves normally in the bulk, but some of the 
mass gets reflected once the mass gets to a boundary. The total mass 
(integral of u over space) then decreases each time it goes through the 
boundary in one direction, and increases if the velocity is set to go the 
other way. Any idea why this might be the case? The periodic boundary works 
fine when its just diffusion and no advection. 

Thank you,
Gary

-- 
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] Obtaining the final position of the vertices

2019-03-19 Thread Wolfgang Bangerth
On 3/19/19 9:57 AM, David F wrote:
> 
> I am not sure how to answer your question. I'm using a very basic a setup 
> equivalent to step-8. Therefore, I have a solution vector with final 
> displacements where each entry corresponds to the displacement of a dof. My 
> aim is to find the initial position of the vertices in the form of 
> std::vector> and the final position (i.e., the initial position + 
> its displacement) in the same form.

I see -- you are using the pre-existing mesh and want to interpret the 
solution as the *displacement* to the mesh.

If all you are interested in is to output the solution on the deformed mesh, 
you can call Data::build_patches() with the mapping argument, where you use a 
MappingQEulerian object initialized with your displacement vector.

If you want to have a representation in the program, you can use the function 
in namespace DoFTools that returns a vector of support point locations for 
each DoF, and then query your solution vector for the corresponding 
displacements. You might have to juggle with the x, y, and z displacements a 
bit to get this right, but the infrastructure is all there.

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] Obtaining the final position of the vertices

2019-03-19 Thread David F
Dear prof. Bangerth, 

I am not sure how to answer your question. I'm using a very basic a setup 
equivalent to step-8. Therefore, I have a solution vector with final 
displacements where each entry corresponds to the displacement of a dof. My 
aim is to find the initial position of the vertices in the form of 
std::vector> and the final position (i.e., the initial position 
+ its displacement) in the same form.


Best,
David.

On Tuesday, 19 March 2019 16:17:58 UTC+1, Wolfgang Bangerth wrote:
>
> On 3/19/19 9:12 AM, David F wrote: 
> > 
> > I want to obtain the final position of the vertices (specifically, the 
> > vertices at the faces), i.e., the deformed configuration. I think that a 
> way 
> > of doing this is by creating a set of points and using fe_values to 
> > extrapolate the solution to those points. However, I don't like the idea 
> of 
> > extrapolating since the displacements are actually known in the vertices 
> > already. Could somebody give me some hint about what's the best way to 
> do this? 
>
> David -- how do you actually deform the mesh? Are you moving around 
> vertices? 
> In that case, you know where you move them to. Or are you using a 
> MappingEulerian? 
>
> 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] Obtaining the final position of the vertices

2019-03-19 Thread Wolfgang Bangerth
On 3/19/19 9:12 AM, David F wrote:
> 
> I want to obtain the final position of the vertices (specifically, the 
> vertices at the faces), i.e., the deformed configuration. I think that a way 
> of doing this is by creating a set of points and using fe_values to 
> extrapolate the solution to those points. However, I don't like the idea of 
> extrapolating since the displacements are actually known in the vertices 
> already. Could somebody give me some hint about what's the best way to do 
> this?

David -- how do you actually deform the mesh? Are you moving around vertices? 
In that case, you know where you move them to. Or are you using a 
MappingEulerian?

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.


[deal.II] Obtaining the final position of the vertices

2019-03-19 Thread David F
Hi all,

I want to obtain the final position of the vertices (specifically, the 
vertices at the faces), i.e., the deformed configuration. I think that a 
way of doing this is by creating a set of points and using fe_values to 
extrapolate the solution to those points. However, I don't like the idea of 
extrapolating since the displacements are actually known in the vertices 
already. Could somebody give me some hint about what's the best way to do 
this?

Thanks in advance.

-- 
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: different results when compute integration with adaptive mesh

2019-03-19 Thread Wolfgang Bangerth
On 2/12/19 12:40 AM, chucui1...@gmail.com wrote:
> 
>  > Constraints are funny and sometimes require deep thought about what exactly
>  > they mean. What happens if you don't apply constraints to the
>  > cell_volume_matrix and cell_gradient_matrix -- i.e., you copy the elements 
> 1:1
>  > into the matrix, without the 'constraints' object? (Like in step-4.)
> I rewrite my code as you say, and the results are:
> |
> Cycle 0:
>     Number of active cells:       20
>     Number of degrees of freedom: 89
>   u square 0: 0.0305621  grad u square 0: 0.256542
> Cycle 1:
>     Number of active cells:       44
>     Number of degrees of freedom: 209
>   u square 1: 0.144303  grad u square 1: 0.49794
> Cycle 2:
>     Number of active cells:       92
>     Number of degrees of freedom: 449
>   u square 2: 0.101366  grad u square 2: 0.477137
> Cycle 3:
>     Number of active cells:       188
>     Number of degrees of freedom: 881
>   u square 3: 0.13305  grad u square 3: 0.529805
> Cycle 4:
>     Number of active cells:       368
>     Number of degrees of freedom: 1737
>   u square 4: 0.145489  grad u square 4: 0.565953
> Cycle 5:
>     Number of active cells:       737
>     Number of degrees of freedom: 3409
>   u square 5: 0.141771  grad u square 5: 0.574863
> Cycle 6:
>     Number of active cells:       1436
>     Number of degrees of freedom: 6705
>   u square 6: 0.171807  grad u square 6: 0.642486
> Cycle 7:
>     Number of active cells:       2729
>     Number of degrees of freedom: 12521
>   u square 7: 0.160819  grad u square 7: 0.625337
> 
> |

Chucui,
you would expect that both the 'u square' and 'grad u square' converge to some 
fixed value, which appears to be the case with the implementation you have here.

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] Computing the curl of a solution vector field obtained from Nedelec elements

2019-03-19 Thread Wolfgang Bangerth
On 2/12/19 6:15 PM, SebG wrote:
> 
> sorry about not getting back to this for while. I have addressed your second 
> point and created an example that is much easier to understand. Now the 
> geometry is the unit cube. The vector field is A = [0 , 0 , y] and then 
> curl(A) = [1 , 0 , 0].
> 
> I have attached three screenshots which compare the analytic result 
> (left-hand 
> side) to the numerical approximation (right-hand side). The approximation of 
> the field, A, seems fine except for a strange wiggle of the magnitude on the 
> finest mesh, see A_magnitude.png. From looking at the graphical output 
> in**curlA_glyphs.png, the approximation of the curl of A does not work well.
> 
> Regarding your third point, I have created the following convergence table:
> 
> cycle cells dofs  field       curl
>    error  rate error    rate
>      0 1    54 0.1866    - 1.4410  -
>      1 8   300 0.1967 0.95 3.0428   0.47
>      2    64  1944 0.1576 1.25 4.8742   0.62
>      3   512 13872 0.1429 1.10 8.8342   0.55
> 
> Regarding the field everything looks good. For the curl, the order seems to 
> be 
> one half although the error increases. This is a bit surprising for me 
> because 
> the overall error increases with the number of dofs. Apart from this, I am 
> asking myself if the errors should not be much smaller because the field is 
> linear and its curl is constant. Hence, these functions should be an element 
> of the approximating finite element space. If I am not wrong, a single 
> element 
> should be sufficient for an approximation in this case. With FE_NedelecSZ the 
> results are similar.

Seb -- as you may have noticed, none of the main developers work on curl 
problems with the Nedelec elements, and so it is difficult to get definitive 
answers about questions.

If you're still chasing answers, let's first focus on the potential (which I 
think the table above labels as "field error". I don't know how the error rate 
is computed, but the numbers for the rates can't be correct: I think they show 
a reduction *factor*, not a rate. The rate is the exponent in
   error = C h^r
where from one line to the next, h is reduced by a factor of 2. So for the 
last two rows, for example, you need to have

   0.1429 = (1/2)^r 0.1576

or r = -log(0.1429/0.1576)/log(2) = 0.14.

So what I see here is that the convergence of the potential (field) is already 
not correct unless you happen to have a very non-smooth right hand side of the 
equation you are trying to solve.

When you later compute the error of the curl of the potential, one typically 
expects that the convergence rate is one lower for the derivative than the one 
for the field itself. This would suggest that you would get a rate of -0.86 -- 
which turns out to be almost exactly what you have in your table above when 
one computes the rate with the formula above.

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: How to compute convergence rate of L2 norm of error without exact solution? And how to compute convergence rate of in different golobally refinement?

2019-03-19 Thread Wolfgang Bangerth
On 2/27/19 2:40 AM, chucui1...@gmail.com wrote:
> 
> I have used SolutionTransfer as you say, but if I set a phi_0 fixed, then 
> project it into finite element space and get a vector phi_0_h, then I get 
> phi_0_h/2, phi_0_h/4, phi_0_h/8 by using SolutionTransfer, but the norm 
> of (phi_0_h- phi_0_h/2), (phi_0_h/2- phi_0_h/4), (phi_0_h/4- phi_0_h/8) are 
> so 
> big, and I cannot get exact value when I compute the convergence rate. So I 
> wander how to deal with it?
> 
> And the attached is the test code, and the results is:
> phi_0_norm_58:  0.335377  phi_0_norm_68:  0.317158  phi_0_norm_78:  0.277344

Chucui,
have you solved this problem since?

It's not quite clear to me what exactly you are doing here, i.e., which of the 
quantities you have above are *computed* on each mesh, and which ones are 
interpolated from the previous mesh.

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] what is the different between VectorTools::point_value and fe_values.get_function_values()?

2019-03-19 Thread Wolfgang Bangerth


> I think that with the H curl conforming element like FeNedelec or 
> FeNedelecSZ, 
> my solution vector has continuous tangential component, isn't it?

Yes.


> Thus, Is there some where I can see my vector solution is discontinuous? Can 
> I 
> check it by visualization?

Yes, you can do that, but you will generally have to use a rather coarse mesh 
and zoom in quite far to make the discontinuity visible to the eye.


> I think my below question is not related to the topic, but may I ask some 
> thing about the shape function of FeNedelec or FeNedelecSZ?

Yes, but why don't you open a new thread for this so that it's easier for 
others to find questions/answers related to their work?

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

2019-03-19 Thread Wolfgang Bangerth
On 3/18/19 4:30 PM, jane@jandj-ltd.com wrote:
> 
> To impose strongly - would you just 
> useVectorTools::compute_nonzero_tangential_flux_constraints with the 
> ZeroFunction?
> or is there a function similar to compute_no_normal_flux_constraints?

Yes, this will compute the constraints that correspond to u.t=0.

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] issues with manifold id (at least in codim-1) revealed by PR 7775

2019-03-19 Thread Wolfgang Bangerth
On 3/19/19 4:17 AM, luca.heltai wrote:
> what would you think of a function in GridTools like
> 
> GridTools::propagate_internal_manifold_ids(tria, disambiguation_function)
> 
> that would loop over all cells, loop over all faces, check if the neighbor 
> manifold id is the same of this cell, and
> i) the two ids are the same: assign the same id to the face (calling 
> cell->face(f)->set_all_manifold_ids())
> ii) the two ids are different: set the id of the face to the return value of 
> the disambiguation_function(id1, id2)
> 
> disambiguation_function could be of type
> 
> std::function )>
> 
> and its default implementation could be
> 
> std::function )> = [](const manifold_id , const manifold_id ){return 
> std::min(id1, id2);}

You'd have to do the same for edges in 3d, where the situation is more 
difficult because there may be an arbitrary number of adjacent cells. So maybe 
the function should take a std::set?

I think that the min of ids is a rather arbitrary choice. Why not just return 
invalid_manifold_id if the two ids are different? (Which could then be 
interpreted as not setting anything for that face/edge -- or maybe that's what 
should be set, I don't know.)

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

2019-03-19 Thread jane . lee
Thanks Jean-Paul, I set n_components to dim and it ran. 

However, no difference whatsoever in the solution to when I was 
equivalently imposing in the weak form (where the tangential term 
disappears due to being zero), so I do continue to wonder whether the two 
are equivalent. Thank you

On Tuesday, March 19, 2019 at 11:35:45 AM UTC, Jean-Paul Pelteret wrote:
>
> Dear Jane,
>
> const Functions::ZeroFunction no_tang_bcs;
>
>
> This should be const Functions::ZeroFunction 
> no_tang_bcs(n_components);
>
> Best,
> Jean-Paul
>
>

-- 
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-03-19 Thread Jean-Paul Pelteret
Dear Jane,

> const Functions::ZeroFunction no_tang_bcs;

This should be const Functions::ZeroFunction no_tang_bcs(n_components);

Best,
Jean-Paul

-- 
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-03-19 Thread jane . lee
In addition to above, I thought I'd try VectorTools::compute_
nonzero_tangential_flux_constraints:

I did:
std::set no_tang_flux_boundaries;
no_tang_flux_boundaries.insert(1);
const Functions::ZeroFunction no_tang_bcs;
typename FunctionMap::type no_tang_map;
no_tang_map[1] = _tang_bcs;
VectorTools::compute_nonzero_tangential_flux_constraints (dof_handler, 0, 
no_tang_flux_boundaries, no_tang_map, constraints);


And this gives me a dimension miss match error - I'm guessing coming from 
the ZeroFunction, somehow, I thought I might need to use a 
fe.component_mask(velocities), but this isn't one of the arguments. 
Where is the dimension miss match coming from?? 



On Monday, March 18, 2019 at 10:30:32 PM UTC, jane...@jandj-ltd.com wrote:
>
> No that's fair enough. 
>
> I had thought the way I was doing it would be the equivalent of setting no 
> tangential stresses. I actually also did it this way as I wasn't sure how 
> you impose it strongly. 
>
> To impose strongly - would you just use 
> VectorTools::compute_nonzero_tangential_flux_constraints 
> with the ZeroFunction?
> or is there a function similar to compute_no_normal_flux_constraints?
>
> I'll try that and see if there's a difference. 
>
> Thanks
>
> On Monday, March 18, 2019 at 8:22:31 PM UTC, Wolfgang Bangerth wrote:
>>
>> On 3/18/19 12:37 PM, jane...@jandj-ltd.com wrote: 
>> > 
>> > For the u_t=0 condition, I had been imposing weak. So basically, I have 
>> > separated a Neumann boundary condition into: 
>> > n.[pI-2e] = (n.[pI-2e]n)n + (n.[pI-2e]t)t 
>> > and saying that the second term on the rhs is 0 so disappears, and the 
>> first 
>> > term is known weakly imposed into 
>> > topstress_values[q_point]*fe_face_values.normal_vector[q_point] etc. 
>> > Is this not a valid way? 
>>
>> I don't know. I don't, because I don't know how you impose the u.t=0 
>> boundary 
>> condition weakly. The way this is typically done is via the Nitsche 
>> method 
>> (which is essentially what DG approaches use for all interior faces as 
>> well). 
>> But I'd have to spend substantially more time than I have thinking 
>> through the 
>> implications for the case you pose here... 
>>
>> 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] issues with manifold id (at least in codim-1) revealed by PR 7775

2019-03-19 Thread luca.heltai
Nicola, 

what would you think of a function in GridTools like 

GridTools::propagate_internal_manifold_ids(tria, disambiguation_function)

that would loop over all cells, loop over all faces, check if the neighbor 
manifold id is the same of this cell, and 
i) the two ids are the same: assign the same id to the face (calling 
cell->face(f)->set_all_manifold_ids())
ii) the two ids are different: set the id of the face to the return value of 
the disambiguation_function(id1, id2)

disambiguation_function could be of type

std::function

and its default implementation could be 

std::function = [](const manifold_id , const manifold_id ){return 
std::min(id1, id2);}

Would this address your needs?

L.


> On 15 Mar 2019, at 15:16, Nicola Giuliani  wrote:
> 
> I agree that the choice is ambiguous on the faces. 
> 
> However a practical problem remains, let's assume that a user wants to import 
> a very big geometry with a lot of internal faces. Right now he would need to 
> specify all the possible faces in the input file (with an increasing 
> possibility of confusing the manifold ids). 
> I agree he should (otherwise is a criminal) specify what happens when two 
> different manifolds meet but if we have a face neighbouring two cells with 
> the same manifold_id couldn't (or shouldn't) we assume the manifold of the 
> face to be the same? 
> Alternatively we could write a utility to transfer the manifold information 
> to the faces, it should be "morally" the same of calling set_all_manifold_ids.

-- 
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] Difference between MeshWorker::mesh_loop and WorkStream::run

2019-03-19 Thread luca.heltai
> Is there a reason then that there are several examples using 
> integration_loop(), but (afaik) only one using mesh_loop?

Yes. mesh_loop is more recent w.r.t. integration_loop.

mesh_loop was introduced to address some of the oddities that are in 
integration_loop, that make its use somewhat less intuitive w.r.t. 
WorkStream::run, and was originally thought as the “kernel” of integration_loop.

The original idea was to replace the internals of integration_loop with 
mesh_loop, but this stage of the project is not yet completed (for lack of time 
of the original developers of mesh_loop… cough cough… :) …).

The examples in 

https://github.com/dealii/dealii/pull/7806

show you a few more ways to use mesh_loop, for DG methods, or using Automatic 
Differentiation.

Hopefully there will be more examples soon on using mesh_loop.

L.


-- 
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] Difference between MeshWorker::mesh_loop and WorkStream::run

2019-03-19 Thread 'Maxi Miller' via deal.II User Group
Is there a reason then that there are several examples using 
integration_loop(), but (afaik) only one using mesh_loop?

Am Montag, 18. März 2019 19:46:41 UTC+1 schrieb Luca Heltai:
>
> Take a look at this PR for a few examples of usage of mesh_loop: 
>
> https://github.com/dealii/dealii/pull/7806 
>
> The function WorkStream::run 
>
> takes (a minimum of) 5 arguments: 
>
> WorkStream::run(cell, endc, cell_worker, copier, scratch, copy); 
>
> initial and final iterator, a worker function, a copier function, a 
> scratch object, and a copy object. 
>
> Logically, WorkStream::run does the following: 
>
> for(it = cell; it != endc; ++it) { 
> cell_worker(cell, scratch, copy); 
> copier(copy); 
> } 
>
> and is agnostic about what type of iterator you pass to it. 
>
> MeshWorker::mesh_loop is a specific implementation of the above in which 
> the iterator is a cell iterator, and the loop is thought to run also on 
> each face of each cell as well, in order to assemble DG type contributions, 
> i.e., it allows you to specify two additional functions, a boundary_worker 
> and a face_worker: 
>
> WorkStream::mesh_loop(cell, endc, cell_worker, copier, scratch, copy, 
> flags, boundary_worker, face_worker); 
>
> This is logically identical to 
>
> for(it = cell; it != endc; ++it) { 
> cell_worker(cell, scratch, copy); 
> for(unsigned int f=0; fif(cell->face(f)->at_boundary() { 
> boundary_worker(cell, f, scratch, copy); 
>} else { 
> // Get subface, neighbor face index, and neighbor subface 
> index 
> const auto [sf, nf, nsf] = get_neighbor_info(cell, f); 
> // Run the face worker 
> face_worker(cell, f, sf, cell->neighbor(f), nf, nsf); 
>} 
> } 
> copier(copy); 
> } 
>
> What mesh_loop does is controlled in a finer way by the assembly flags. 
> For example it is possible to visit each face only once, or to assemble 
> first faces and then cells, etc. 
>
> MeshWorker::integration_loop 
>
> is a specific implementation of MeshWorker::mesh_loop where you use the 
> DoFInfo and DoFInfoBox objects, that contain scratch and copy data with a 
> different structure. The two behave in a substantially identical way, but 
> with different data structures. 
>
> In particular, you could implement MeshWorker::integration_loop using 
> MeshWorker::mesh_loop, and some wrapping around DoFInfo and DoFInfoBox 
> objects. 
>
> I hope this clarifies a little the scope of run (agnostic about iterators) 
> and mesh_loop (actually expecting cell iterators, and taking care of 
> subface/face/neighbor matching when looping over faces). 
>
> Best, 
> Luca. 
>
>
> > On 18 Mar 2019, at 17:46, 'Maxi Miller' via deal.II User Group <
> dea...@googlegroups.com > wrote: 
> > 
> > Similar question: How are the different loop-functions differentiated, 
> i.e. MeshWorker::integration_loop and MeshWorker::mesh_loop? Both are able 
> to loop over faces, boundaries and cells, but what are the differences 
> here? 
> > Thanks! 
> > 
> > Am Dienstag, 22. Januar 2019 16:56:28 UTC+1 schrieb Bruno Turcksin: 
> > Le mar. 22 janv. 2019 à 10:48, 'Maxi Miller' via deal.II User Group 
> >  a écrit : 
> > > I. e. if I would like to calculate f. ex. the L2-norm of a vector 
> (while neglecting that there already is a function for that), I can use 
> WorkStream for parallelization of that, but not MeshWorker, is that 
> correct? 
> > That's right. You can use WorkStream using the iterators for the 
> > vector but MeshWorker won't work because it expects cell iterators. 
> > 
> > Best, 
> > 
> > Bruno 
> > 
> > -- 
> > 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+un...@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.