Re: [deal.II] what is the different between VectorTools::point_value and fe_values.get_function_values()?

[Phạm Ngọc Kiên]

Thank you very much for your guidance.
I think that with the H curl conforming element like FeNedelec or
FeNedelecSZ, my solution vector has continuous tangential component, isn't
it?
Thus, Is there some where I can see my vector solution is discontinuous?
Can I check it by visualization?

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?
My testing codes showed that with those elements, the shape functions
changed the values when the mesh grid changed.
For example, assuming I have a cell that is a cube, whose edge length is 1,
the shape function at the edge is {0, 0, 1}.
For the smaller cell whose edge length is 0.5, the shape function at the
edge is {0, 0, 2}.
I saw it was double every time I reduce the edge length by the factor of 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 my solution change when refining the mesh.
I think there should be some idea that we change the shape function values
like this.
Could you please give me some advice?

Best regards,
Pham Ngoc Kien


Vào Th 3, 19 thg 3, 2019 vào lúc 02:22 Wolfgang Bangerth <
bange...@colostate.edu> đã viết:

> On 3/18/19 4:14 AM, Daniel Arndt wrote:
> >
> > VectorTools::point_value() does more or less the same you are doing
> above
> > manually, i.e. calling GridTools::find_active_cell_around_point() and
> then
> > initializing a FEValues object for evaluating the given finite element
> vector.
>
> Specifically, you can take a look here:
>
> https://github.com/dealii/dealii/blob/master/include/deal.II/numerics/vector_tools.templates.h#L8633
>
> 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.
>

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

[jane . lee]

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

[Wolfgang Bangerth]

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: 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] deal.II Newsletter #71

[Rene Gassmoeller]

Hello everyone!

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

#7825: Add missing #ifdef DEAL_II_WITH_THREADS around tbb calls (proposed by 
Rombur) https://github.com/dealii/dealii/pull/7825

#7822: Perform p-refinement during 'execute_coarsening_and_refinement()'. 
(proposed by marcfehling) https://github.com/dealii/dealii/pull/7822

#7820: Check SCALAPACK symbols (proposed by masterleinad) 
https://github.com/dealii/dealii/pull/7820

#7819:  DoFRenumbering::cell_wise with p::Tria (proposed by davydden) 
https://github.com/dealii/dealii/pull/7819

#7818: Set cell manifold ids when the apply_all_indicators_to_manifold flag is 
used. (proposed by masterleinad) https://github.com/dealii/dealii/pull/7818

#7817: hp::DoFHandler::execute_coarsening_and_refinement() (proposed by 
marcfehling) https://github.com/dealii/dealii/pull/7817

#7816: CMake: Reorganize MPI sanity check (proposed by tamiko) 
https://github.com/dealii/dealii/pull/7816

#7815: Only permit use of ld.lld linker when the Clang compiler is used. 
(proposed by jppelteret; merged) https://github.com/dealii/dealii/pull/7815

#7814: Fix 'test' target for MSVC (proposed by masterleinad) 
https://github.com/dealii/dealii/pull/7814

#7813: Add some basic tests for SymEngine (proposed by jppelteret) 
https://github.com/dealii/dealii/pull/7813

#7810: Add CMake module and configuration files for SymEngine (proposed by 
jppelteret; merged) https://github.com/dealii/dealii/pull/7810

#7808: Encapsulate the calls to tbb::parallel_for (proposed by Rombur; merged) 
https://github.com/dealii/dealii/pull/7808

#7807: Fix compiling bundled boost with MSVC and C++17 (proposed by 
masterleinad) https://github.com/dealii/dealii/pull/7807

#7806: "Reasonable" ScratchData and CopyData (proposed by luca-heltai) 
https://github.com/dealii/dealii/pull/7806

#7805: Improve the error text of an assertion. (proposed by bangerth; merged) 
https://github.com/dealii/dealii/pull/7805

#7804: remove unneeded assert in internal::update_ghost_values_finish() 
(proposed by davydden; merged) https://github.com/dealii/dealii/pull/7804

#7803: Fix MatrixFree::loop() for the case where VectorType is not a vector 
(proposed by kronbichler; merged) https://github.com/dealii/dealii/pull/7803

#7801: Allow CXX14/17 autodetection for MSVC (proposed by masterleinad) 
https://github.com/dealii/dealii/pull/7801

#7800: Update allowed MSVC version in bundled boost (proposed by masterleinad; 
merged) https://github.com/dealii/dealii/pull/7800

#7799: Update documentation of class Vector. (proposed by bangerth; merged) 
https://github.com/dealii/dealii/pull/7799

#7798: Remove make slice (proposed by drwells; merged) 
https://github.com/dealii/dealii/pull/7798

#7797: Fix warnings in Convert::to_pattern() (proposed by masterleinad; merged) 
https://github.com/dealii/dealii/pull/7797

#7795: enable lld linker if available (proposed by tjhei; merged) 
https://github.com/dealii/dealii/pull/7795

#7793: Add useful methods to TrilinosWrappers::PreconditionAMG::AdditionalData 
(proposed by jppelteret) https://github.com/dealii/dealii/pull/7793

#7792: Fix link targets in frames. (proposed by drwells) 
https://github.com/dealii/dealii/pull/7792

#7791: Fix renumbering of level dofs for FE_Q (proposed by kronbichler; merged) 
https://github.com/dealii/dealii/pull/7791

#7787: add and use type traits for data exchange in matrix-free (proposed by 
davydden; merged) https://github.com/dealii/dealii/pull/7787

#7785: Fix a bug in cuda add_and_dot (proposed by Rombur; merged) 
https://github.com/dealii/dealii/pull/7785

#: Introduced p-refinement and p-coarsening flags. (proposed by 
marcfehling; merged) https://github.com/dealii/dealii/pull/

#7776: Signals for p::d::Triangulation::repartition. (proposed by marcfehling; 
merged) https://github.com/dealii/dealii/pull/7776

#7760: Fix compiling with Sacado and CUDA (proposed by masterleinad; merged) 
https://github.com/dealii/dealii/pull/7760

#7751: interpolate_to_different_mesh for hp::DoFHandler (proposed by 
starki0815; merged) https://github.com/dealii/dealii/pull/7751


## And this is a list of recently opened or closed discussions:

#7824: Undeclared identifier 'tbb' (opened) 
https://github.com/dealii/dealii/issues/7824

#7823: Document CUDA support better, write a tutorial step (opened) 
https://github.com/dealii/dealii/issues/7823

#7821: AMD GPU support (opened) https://github.com/dealii/dealii/issues/7821

#7812: Edit step-22 InverseMatrix discussion (opened) 
https://github.com/dealii/dealii/issues/7812

#7811: ld.lld cannot link with openmpi (with enabled fortran support in 
openmpi) (opened) https://github.com/dealii/dealii/issues/7811

#7809: deal.II Release 9.1 5/17/2019 (opened) 
https://github.com/dealii/dealii/issues/7809

#7802: read_ucd ignores manifold ids 

Re: [deal.II] Difference between MeshWorker::mesh_loop and WorkStream::run

[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; fface(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 
>  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+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.

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

[jane . lee]

Hi Wolfgang,

step-20: Yes indeed I do agree that that is what I am doing. I guess I'm 
trying now to find out what else it could be that is producing: the correct 
boundary points as in the Dirichlet condition when local_rhs=0 is done 
again (overwriting), but the wrong boundary points when it isn't. 


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?
(PS apologies - earlier, i had said my boundary condition is partial but 
i've imposed it in a Neumann sense. sorry for any confusion)


On Monday, March 18, 2019 at 5:17:08 PM UTC, Wolfgang Bangerth wrote:
>
>
> Jane, 
>
> > Ok, re the first step-20 bc issue, I'll have another think, but I still 
> am not 
> > sure why then it isn't giving me the exact figure, whilst my suggestion 
> is (I 
> > understand your point here and I would have said I agreed with you, but 
> my 
> > implementation does work). I've verified this part of the code using 
> exact 
> > solutions so I really have no idea what the issue is, but I guess i'll 
> have to 
> > have another think. 
>
> In your exact solution, you may simply have had zero boundary conditions 
> on 
> parts of the boundary, so removing these contributions again made no 
> difference. 
>
> I don't know if that rings any kind of bell, but if you're in 2d and you 
> have 
> a rectangular geometry, you first visit the left and right boundaries of 
> cells, and then the bottom and top boundaries. So if you do local_rhs=0 
> before 
> every boundary face, in essence you are overwriting everything you've done 
> for 
> the left/right boundaries if you are on a cell that also has a bottom or 
> top 
> boundary. 
>
>
> > On the second, step-22 issue, i don't have u.n=0 conditions. I have 
> partial 
> > boundary conditions in the form of that stated in the tutorial (1st in 
> the 
> > fourth set of Bc explanations) of: 
> > u_t = 0 
> > n.(n.(pI-2*epsilon)) = nonzero value. 
> > So this is imposed weakly, and no strong conditions are set on that 
> boundary 
> > at all, and that nonzero value equals the topstress_value I had before 
> in the 
> > weak form. 
>
> But u_t=0 is a strong boundary condition. Are you not calling a function 
> that 
> imposes it? 
>
> 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] what is the different between VectorTools::point_value and fe_values.get_function_values()?

[Wolfgang Bangerth]

On 3/18/19 4:14 AM, Daniel Arndt wrote:
> 
> VectorTools::point_value() does more or less the same you are doing above 
> manually, i.e. calling GridTools::find_active_cell_around_point() and then 
> initializing a FEValues object for evaluating the given finite element 
> vector. 

Specifically, you can take a look here:
https://github.com/dealii/dealii/blob/master/include/deal.II/numerics/vector_tools.templates.h#L8633

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] projecting boundary values for H(div) [and maybe H(curl)]

[Wolfgang Bangerth]

On 3/16/19 7:12 AM, Konrad wrote:
> 
> Tried it with BDM elements but it does not work. I strongly think about 
> starting to contribute to deal.ii.

I'm glad to see that you've already found the issue, but I would like to add 
that I think this here would be a really good idea! :-) We are always happy to 
have new contributors!

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

[Wolfgang Bangerth]

Jane,

> Ok, re the first step-20 bc issue, I'll have another think, but I still am 
> not 
> sure why then it isn't giving me the exact figure, whilst my suggestion is (I 
> understand your point here and I would have said I agreed with you, but my 
> implementation does work). I've verified this part of the code using exact 
> solutions so I really have no idea what the issue is, but I guess i'll have 
> to 
> have another think.

In your exact solution, you may simply have had zero boundary conditions on 
parts of the boundary, so removing these contributions again made no difference.

I don't know if that rings any kind of bell, but if you're in 2d and you have 
a rectangular geometry, you first visit the left and right boundaries of 
cells, and then the bottom and top boundaries. So if you do local_rhs=0 before 
every boundary face, in essence you are overwriting everything you've done for 
the left/right boundaries if you are on a cell that also has a bottom or top 
boundary.


> On the second, step-22 issue, i don't have u.n=0 conditions. I have partial 
> boundary conditions in the form of that stated in the tutorial (1st in the 
> fourth set of Bc explanations) of:
> u_t = 0
> n.(n.(pI-2*epsilon)) = nonzero value.
> So this is imposed weakly, and no strong conditions are set on that boundary 
> at all, and that nonzero value equals the topstress_value I had before in the 
> weak form.

But u_t=0 is a strong boundary condition. Are you not calling a function that 
imposes it?

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

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

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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] projecting boundary values for H(div) [and maybe H(curl)]

[Konrad]

Thank you, Daniel. That confirms my intuition about that. The push forward 
from the ref cell to physical cells is, I found after thinking about it 
again, for RaviartThomas, BDM, Nedelec elements etc not the same as a Piola 
transform (which eliminates scaling but not rotation). My code works now 
without having to rescale.

Thanks again :-)

Konrad

On Monday, March 18, 2019 at 11:17:40 AM UTC+1, Daniel Arndt wrote:
>
>  Konrad,
>  
>  I implemented a class derived from the Function class that evaluates 
> a scalar or vector shape function at a given (set of) point(s) in a 
> physical cell. If I check the output graphically I see that the vector 
> shapefunctions on the physical cell, for example for lowest order 
> Raviart-Thomas elements has a magnitude that scales with the inverse edge 
> length. But on any physical cell the magnitude should be one I believe.
>  
>   For the Raviart-Thomas the degrees of freedom for the lowest order 
> element should be the face integrals. Hence, I would expect them to scale 
> with the (inverse of the) size of the face.
>  
>  
> I assume that has something to do with the mapping used by the FEValues 
> class to map to the physical cell, right? Does the FEValues class initiated 
> with any vector element use the Piola transform? If not how can I influence 
> this?
>  
>  Have a look at 
> https://www.dealii.org/current/doxygen/deal.II/group__mapping.html. The 
> transformation type used in Mapping::transform() is normally set by the 
> given FiniteElement class.
>  
>  Best,
>  Daniel
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Finding the y-coordinate of the boundary during the assembly

[jane . lee]

This did exactly what I needed. Thanks!!

On Saturday, March 16, 2019 at 5:48:14 PM UTC, Jean-Paul Pelteret wrote:
>
> There’s a function to compute the bounding box in GridTools:
>
> https://www.dealii.org/current/doxygen/deal.II/namespaceGridTools.html#ae1ec55abefa31cf001fd29d8d4d993f1
>
>
> On 16 Mar 2019, at 18:36, jane...@jandj-ltd.com  wrote:
>
> Thanks for your suggestions.
>
> BoundingBox should for now be sufficient, but I am having some trouble 
> using it properly. 
>
> I am unsure what to put in to use it? triangulation? cell? 
>
> I think i need something like
> BoundingBox mybox;
> mybox.get_boundary_points();
>
> but I'm a little unsure where I put my argument and what that is - I can't 
> seem to find any examples either. Sorry if this is very basic!
>
> On Friday, March 15, 2019 at 11:40:28 PM UTC, Wolfgang Bangerth wrote:
>>
>>
>> > Yes, what you put into much better words than mine is exactly what I am 
>> > needing - For a given quadrature point at 
>> > (x,y), find how far the domain extends above (x,y) in y-direction? 
>> > 
>> > So I am looking to find the y-coordinate of the point which is directly 
>> > above the (x,y) in question, so that I can find how far it extends 
>> above 
>> > it (by subtracting it from what I am trying to find) 
>>
>> As Jean-Paul already mentioned, this is easy enough if your top surface 
>> is level (i.e., at the same y value). It is, in general, difficult to 
>> figure this out on unstructured meshes if the top surface is not level. 
>>
>> The way to do this then is to define a "depth" variable D(x,y). You know 
>> that D(x,y)=0 at the top surface, and that it grows linearly with depth 
>> y. So D(x,y) satisfies a differential equation of the form 
>>
>>d/d(-y) D(x,y) = 1 
>>
>> or equivalently 
>>
>>d/dy D(x,y) = -1 
>>
>> which you can also write as follows: 
>>
>>   (0,-1) . nabla D(x,y) = 1 
>>
>> So it is an advection equation with the advection velocity being the 
>> vector pointing straight down. The boundary condition at the "inflow" 
>> boundary -- which here is the top boundary) is D(x,y)=0. 
>>
>> Then, if you need to know the depth at a given point (x_q,y_q), for 
>> example at a quadrature point when assembling your linear system, all 
>> you need to do is evaluate the depth field D(x_q,y_q) that you have 
>> previously computed. 
>>
>> 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+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.

Re: [deal.II] projecting boundary values for H(div) [and maybe H(curl)]

[Daniel Arndt]

 Konrad,
 
 I implemented a class derived from the Function class that evaluates 
a scalar or vector shape function at a given (set of) point(s) in a 
physical cell. If I check the output graphically I see that the vector 
shapefunctions on the physical cell, for example for lowest order 
Raviart-Thomas elements has a magnitude that scales with the inverse edge 
length. But on any physical cell the magnitude should be one I believe.
 
  For the Raviart-Thomas the degrees of freedom for the lowest order 
element should be the face integrals. Hence, I would expect them to scale 
with the (inverse of the) size of the face.
 
 
I assume that has something to do with the mapping used by the FEValues 
class to map to the physical cell, right? Does the FEValues class initiated 
with any vector element use the Piola transform? If not how can I influence 
this?
 
 Have a look at 
https://www.dealii.org/current/doxygen/deal.II/group__mapping.html. The 
transformation type used in Mapping::transform() is normally set by the 
given FiniteElement class.
 
 Best,
 Daniel

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] what is the different between VectorTools::point_value and fe_values.get_function_values()?

[Daniel Arndt]

>
> [...] 
>
1.  When evaluating point is on the vertex of a cell, this means that there 
> should be several cells in 3D have this vertex. 
> Thus, my first question is the function :
>
> const std::pair::active_cell_iterator, Point>
> cell_point = GridTools::find_active_cell_around_point(mapping, 
> dof_handler, point); 
>
> I personally think that as I set it to be const, so the return cell 
> (cell_point.first) is only one cell, rather than a group of cells. Could you 
> please tell me about this?
>
>
This function gives you one cell that includes the given point, not all of 
them. In case, you are using continuous elements, the value you obtain 
should not depend on the choice of this cell. Otherwise, you have to be 
more specific which cell you want.
 

> 2. My second idea is with the codes:
> const Quadrature 
> quadrature(GeometryInfo::project_to_unit_cell(cell_point.second));
> FEValues fe_values(mapping, dof_handler.get_fe(), quadrature,
> update_values | update_gradients | 
> update_quadrature_points);
> fe_values.reinit(cell_point.first);
>
> I have only one quadrature defined on the cell, and the printing result show 
> me the same thing. 
> Can I use this way to get the solution, or it would be more efficient when 
> using codes like:
> QGaussLobatto quadrature_formula(2);
> FEValues fe_values(fe, quadrature_formula,
> update_values | update_gradients |
> update_quadrature_points | update_JxW_values);
> fe_values.reinit(cell_point.first);
>
> And then finding the quadrature point coincide with the evaluating point to 
> get the solution?
>
>
If you know beforehand that your evaluation coincides with a unit support 
point on a known cell, you can save the call to 
GridTools::find_active_cell_around_point().
Otherwise, you likely can't do much better than using its return value.
 

> 3. My third thought is in my function I have only 1 quadrature point so that 
> temp_solution_re[0] and temp_solution_im[0] are the variables containing the 
> solution. 
>
> The printing results showed that the solution from point_value() is identical 
> with those from temp_solution_re[0] and temp_solution_im[0].
>
> Do temp_solution_re[1], temp_solution_re[2], 
> temp_solution_im[1],temp_solution_im[2] make any sense for my solution? 
> Because I saw the code can run with these variables.
>
>
VectorTools::point_value() does more or less the same you are doing above 
manually, i.e. calling GridTools::find_active_cell_around_point() and then 
initializing a FEValues object for evaluating the given finite element 
vector. Hence, it is not surprising if the results coincide. Since the size 
of temp_solution_re and temp_solution_im is 1, it doesn't make sense to 
access the values in the second and third entry.

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: issues with manifold id (at least in codim-1) revealed by PR 7775

[Daniel Arndt]

Nicola,

-) I believe that GridIn.read_ucd appears to ignore manifold ids associated 
> to the cells, it seems to me that the flag 
> apply_all_indicators_to_manifold_ids does not affect the manifold id for 
> the cell. This issue has never been seen because internally the material_id 
> was used.
>

Can you confirm that https://github.com/dealii/dealii/pull/7818 fixes this 
problem for you?

-) As far as I understand the default behaviour is to consider only user 
> specified manifold ids for the lines of a cell otherwise flat_manifold_id. 
> Would it make sense to use the cell->manifold_id instead of 
> flat_manifold_id?
>

Note that you can still specify the material id/manifold id for faces and 
lines in the input grid file.

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] New FiniteElement; Cohesive-zone model; A-FEM

[luca.heltai]

You may also find useful this link:

https://journals.ub.uni-heidelberg.de/index.php/ans/issue/view/2930

where the authors implement XFEM in deal.II to model interface problems.

L.

> On 18 Mar 2019, at 9:39, Daniel Arndt  
> wrote:
> 
> James,
>  
> I am a graduate student who will be defending in the not-too-distant future, 
> and afterwards I want to remain in the research community even though I will 
> be leaving academia for the foreseeable future. My research group has an 
> established cohesive-zone model (CZM) code that is used in Abaqus, and I have 
> been using this for all of my graduate work. The main drawback to not 
> switching over to deal.ii years ago was that it lacks the ability to model 
> crack-propagation (which is the main focus of my research), and the 
> resistance to changing from a code the group knows already works.
>  
> There are at least some papers that use deal.II for crack propagations, e.g. 
> "A primal-dual active set method and predictor-corrector mesh adaptivity for 
> computing fracture propagation using a phase-field approach".
> 
> 
> With this in mind, would a user/developer be willing to guide me along the 
> path of writing a new FiniteElement? I have already written initial 
> traction-separation laws (in C++) that could be used in a "volume" 
> cohesive-zone element. I have collected papers of others who have already 
> written cohesive-zone finite-elements for the commercial software package 
> Abaqus, and my intention would be to replicate this work (with perhaps some 
> small modifications/improvements) for deal.ii so that I can continue work in 
> fracture mechanics.
>  
> It is probably best if you describe what you need and how your new finite 
> element would look like. In general, you need to define a suitable 
> (polynomial) ansatz space and how the degrees of freedom are evaluated. In 
> case, the element uses tensor-product polynomials there are alreay at lot of 
> prerequsites available. You might one want to look into the implementation of 
> FE_Q (https://github.com/dealii/dealii/blob/master/source/fe/fe_q.cc)
> 
>  
> Additionally, does deal.ii have the capabilities of the so-called A-FEM 
> (augmented finite element method)? A-FEM allows one element to change an 
> element from one finite-element into another finite-element during the 
> calculation. An example from fracture mechanics: starting with an initially 
> isotropic elastic calculation of generic shape, when a certain stress 
> threshold is met by an element, the elastic element is replaced with a CZM 
> element to allow for arbitrary crack-initiation and propagation to occur. 
> Note: other FEM-based methods appear to use this same idea of replacing one 
> element type with another during a calculation, but I am do not have 
> extensive knowledge of the nuanced differences between them.
>  
> deal.II support hp-adaptivity that requires the possibility to use different 
> finite elements across the discretization. Have a look at 
> https://www.dealii.org/current/doxygen/deal.II/step_27.html.
> 
> Best,
> Daniel
> 
> -- 
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see 
> https://groups.google.com/d/forum/dealii?hl=en
> --- 
> You received this message because you are subscribed to the Google Groups 
> "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an 
> email to dealii+unsubscr...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: New FiniteElement; Cohesive-zone model; A-FEM

[Daniel Arndt]

James,
 

> I am a graduate student who will be defending in the not-too-distant 
> future, and afterwards I want to remain in the research community even 
> though I will be leaving academia for the foreseeable future. My research 
> group has an established cohesive-zone model (CZM) code that is used in 
> Abaqus, and I have been using this for all of my graduate work. The main 
> drawback to not switching over to deal.ii years ago was that it lacks the 
> ability to model crack-propagation (which is the main focus of my 
> research), and the resistance to changing from a code the group knows 
> already works.
>
 
There are at least some papers that use deal.II for crack propagations, 
e.g. "A primal-dual active set method and predictor-corrector mesh 
adaptivity for computing fracture propagation using a phase-field approach 
".


With this in mind, would a user/developer be willing to guide me along the 
> path of writing a new FiniteElement? I have already written initial 
> traction-separation laws (in C++) that could be used in a "volume" 
> cohesive-zone element. I have collected papers of others who have already 
> written cohesive-zone finite-elements for the commercial software package 
> Abaqus, and my intention would be to replicate this work (with perhaps some 
> small modifications/improvements) for deal.ii so that I can continue work 
> in fracture mechanics.
>
 
It is probably best if you describe what you need and how your new finite 
element would look like. In general, you need to define a suitable 
(polynomial) ansatz space and how the degrees of freedom are evaluated. In 
case, the element uses tensor-product polynomials there are alreay at lot 
of prerequsites available. You might one want to look into the 
implementation of FE_Q 
(https://github.com/dealii/dealii/blob/master/source/fe/fe_q.cc)

 

> Additionally, does deal.ii have the capabilities of the so-called A-FEM 
> (augmented finite element method)? A-FEM allows one element to change an 
> element from one finite-element into another finite-element during the 
> calculation. An example from fracture mechanics: starting with an initially 
> isotropic elastic calculation of generic shape, when a certain stress 
> threshold is met by an element, the elastic element is replaced with a CZM 
> element to allow for arbitrary crack-initiation and propagation to occur. 
> Note: other FEM-based methods appear to use this same idea of replacing one 
> element type with another during a calculation, but I am do not have 
> extensive knowledge of the nuanced differences between them.
>
 
deal.II support hp-adaptivity that requires the possibility to use 
different finite elements across the discretization. Have a look at 
https://www.dealii.org/current/doxygen/deal.II/step_27.html.

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.