Re: [deal.II] point_value FE_Nedelec

2021-04-18 Thread John Smith
 

Yes, I have. Thus the question.

I would like to have a simple silly test to see that the FE_Nedelec 
elements are what I expect them to be. So, I made a simple program similar 
to step-3 that sets up a test. The program uses FE_Nedelec<2>. It does all 
the steps similarly to step-3. However, it does not solve the system of 
linear equations. It substitutes all zeros into the solution vector and 
then makes solution(0) = 1.

Then I have fed this solution together with the dof handler to the 
VectorTools::point_value on a grid of regularly spaced points. The results 
were written into a *.csv file. Then I have made the vector plot out of it. 
The 0-th order, FE_Nedelec<2> fe(0), looks just as what I would expect. The 
2-nd order, however, looks just like the 0-th order. I opened the files in 
spreadsheet and indeed – the shape functions of the 0-th and 2-nd order are 
the same. I ruled out the possibility of mistake (the program prints 
fe.degree) and began worrying. I think VectorTools::point_value is to blame 
probably… 

Is there a replacement for it? 

Best, 

John


On Monday, April 19, 2021 at 4:59:25 AM UTC+2 Wolfgang Bangerth wrote:

> On 4/18/21 7:43 PM, John Smith wrote:
> > 
> > Is VectorTools::point_value supposed to work with FE_Nedelec?
> > 
> > The documentation looks very promising:
> > 
> > Evaluate a possibly vector-valued finite element function defined by the 
> given 
> > DoFHandler 
> > <
> https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.dealii.org%2Fcurrent%2Fdoxygen%2Fdeal.II%2FclassDoFHandler.html=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7C39054f37290545641e5608d902d486d0%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637543934089426552%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000=cfNPiGUOk%2BsIlssgAev6wBAmMiD%2BsZHYrl32blPab0E%3D=0>
>  
>
> > and nodal vector fe_function at the given point point, and return the 
> (vector) 
> > value of this function through the last argument.
>
> John -- have you tried? Oftentimes, just giving it a try is going to 
> answer 
> your question :-)
> Best
> W.
>
>
> -- 
> 
> Wolfgang Bangerth email: bang...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/3cfc8137-d48e-4b98-8870-3dab4f64365en%40googlegroups.com.


Re: [deal.II] Re: interpolate FE_Nelelec

2021-04-18 Thread Wolfgang Bangerth

On 4/18/21 12:24 PM, John Smith wrote:


Thank you for your reply. It is a bit difficult to read formulas in text. So I 
have put a few questions I have in a pdf file. Formulas are better there. It 
is attached to this message. May I ask you to have a look at it?


John:

If I understand your first question right, then you are given a vector field J 
and you are looking for a vector field T so that

  curl T = J
I don't really have anything to offer to this kind of question. There are many 
vector fields T that satisfy this because of the large null space of curl. You 
have a number of condition you would like to "approximate" but there are many 
ways to "approximate" something. In essence, you have two goals: To satisfy 
the equation above and to approximate something. You have to come up with a 
way to weigh these two goals. For example, you could look to minimize

  F(T) = ||curl T - J||^2 + alpha ||T-T_desired||^2
where T_desired is the right hand side of (0.0.5) and you have to determine 
what alpha should be.


As for your other questions: (0.0.9) is indeed called "interpolation"

The issue with the Nedelec element (as opposed to a Q(k)^dim) field is that 
for the Nedelec element,

  \vec phi(x_j)
is not independent of
  \vec phi(x_k)
and so you can't choose the matrix in (0.0.8) as the identity matrix. You 
realize this in (0.0.13). The point I was trying to make is that (0.0.13) 
cannot be exactly satisfied if T(x_j) is not a vector parallel to \hat e_j, 
which in general it will not be. You have to come up with a different notion 
of what it means to solve (0.0.13) because you cannot expect the left and 
right hand side to be equal.


Best
 W>

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/751ebf54-1dc5-fc85-6a99-d898a86b0857%40colostate.edu.


Re: [deal.II] point_value FE_Nedelec

2021-04-18 Thread Wolfgang Bangerth

On 4/18/21 7:43 PM, John Smith wrote:


Is VectorTools::point_value supposed to work with FE_Nedelec?

The documentation looks very promising:

Evaluate a possibly vector-valued finite element function defined by the given 
DoFHandler 
 
and nodal vector fe_function at the given point point, and return the (vector) 
value of this function through the last argument.


John -- have you tried? Oftentimes, just giving it a try is going to answer 
your question :-)

Best
 W.


--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/020523f6-dbe5-82b0-7648-2f6f72348190%40colostate.edu.


[deal.II] point_value FE_Nedelec

2021-04-18 Thread John Smith
 

Hello,

Is VectorTools::point_value supposed to work with FE_Nedelec?

The documentation looks very promising:

Evaluate a possibly vector-valued finite element function defined by the 
given DoFHandler 
 and 
nodal vector fe_function at the given point point, and return the (vector) 
value of this function through the last argument.

Best, 

John.

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/ec5bb84c-e3cb-4f91-91d5-eac424d68d87n%40googlegroups.com.


[deal.II] Re: deal.II-based FSI matrix-free solver / looking for a post-doc

2021-04-18 Thread blais...@gmail.com
Dear Michal,
We have a post-doc position available in my group regarding deal.II and 
topology optimisation for conformal cooling.
Please send an email to bruno.bl...@polymtl.ca if you are interested :)!
Best
Bruno


On Saturday, April 17, 2021 at 6:44:44 p.m. UTC-4 mtwich...@gmail.com wrote:

> Dear all,
>  I have developed a parallel matrix-free monolithic FSI solver in the ALE 
> frame of reference. Exemplary results [1 , 2 
> ], more details in my PhD thesis [3] 
> 
> .
>
> *Highlights*:
>
> -> The proposed monolithic predictor-corrector time integration scheme 
> requires solving the generalized Stokes problem with discontinuous variable 
> coefficients. This has to be done one or two times per time step, depending 
> on the variant of the scheme.
>
> -> I have introduced a new multilevel preconditioner, for a 3D case the 
> time required to solve a system consisting of 421 million unknowns using 96 
> cores is 178 seconds.
>
> -> As demonstrated by the movie[1], I am able to solve 3D problems at 
> moderate Reynolds numbers (time step size of 0.01 did not result in any 
> stability issues with Re=2k). 
>
> ->The preconditioner from the thesis can still be greatly improved. I have 
> a proof of concept demonstrating the basic properties of a new idea of 
> smoother. The results look very promising and I'm curious how it would 
> work. Moreover, the algorithm is suitable for GPU that could result in 
> further speedup.
>  
> [1] 3D "Turek benchmark" at Re = 2000 (30M DoF)
> https://youtu.be/7TaC1dTfAQY
> [2] Water hammer in elastic tube (12M DoF)
> https://youtu.be/TRnlt_eC32Q
>
> Those results were obtained as a part of my PhD thesis. I have finally 
> managed to submit it in January, the defence will be held in 3 weeks 
> (expect reviews next week).
> I am looking for a post-doc (preferably deal.II-related). I'll be happy to 
> get an opportunity to continue my work.
>
> Best,
> Michał Wichrowski
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/e9878b7b-4191-4f1c-9a84-b4cbbb307a6an%40googlegroups.com.


Re: [deal.II] Re: interpolate FE_Nelelec

2021-04-18 Thread John Smith
Dear Wolfgang,

Thank you for your reply. It is a bit difficult to read formulas in text. 
So I have put a few questions I have in a pdf file. Formulas are better 
there. It is attached to this message. May I ask you to have a look at it? 

Best,
John

On Friday, April 16, 2021 at 10:26:51 PM UTC+2 Wolfgang Bangerth wrote:

> On 4/16/21 11:36 AM, John Smith wrote:
> > However, the FE_Nedelec is different. It implements edge elements. They 
> are 
> > vector-based. That is, functions are represented by a superposition of 
> > vector-valued shape functions:
> > 
> > \vec{A} = \sum u_i vec{N}_i .
> > 
> > Therefore, the output of the "value" method in “class CustomFunction : 
> public 
> > Function” must be vector-valued. Three components in a three-dimensional 
> > space.  Otherwise, there is no point in interpolation.
> > 
> > This kind of vectorial approximations is a bread-and-butter topic in 
> > magnetics. See, for example, equation (7) in:
> > 
> > https://ieeexplore.ieee.org/document/497322 
> > <
> https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fieeexplore.ieee.org%2Fdocument%2F497322=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cf7f578c8c6fc44890c3308d900fe30ef%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C1%7C637541914012798794%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000=tJVhiTO272UFKlS0mWzTCrVp1VzPV3zbr2I16KxHkFI%3D=0
> >
> > 
> > In short, the source, i.e the current vector potential, must be 
> projected on 
> > the space spanned by the vector-valued shape functions. Otherwise, the 
> > simulation is numerically unstable. The last, from my experience, is 
> > definitely true.
>
> I think you already found your solution, but just for clarity: When we use 
> the 
> term "interpolation", we typically ask for (scalar) coefficients U_i so 
> that
>
> u_h(x_j) = sum U_i \phi_i(x_j) = g(x_j)
>
> where g is given and x_j are the node points. The problem is that for the 
> Nedelec element, this is not always possible: Not all possible vectors 
> g(x_j) 
> can be represented. This makes sense because if you have N node points in 
> 3d, 
> then you have N scalar coefficients U_i, but you have 3N components of the 
> values g(x_j).
>
> One way to deal with this is to associate a vector, let's say y_j with 
> every 
> node point x_j, and require that
>
> y_j \cdot u_h(x_j) = sum U_i y_j \cdot \phi_i(x_j) = y_j \cdot g(x_j)
>
> These y_j could, for example, be the tangential direction associated with 
> the 
> shape function phi_j. One *could* call this a variation of the term 
> "interpolation", but it is not what VectorTools::interpolate() implements. 
> It 
> would probably not be terribly difficult to implement this kind of 
> function, 
> however!
>
> Best
> W.
>
>
>
> -- 
> 
> Wolfgang Bangerth email: bang...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/4149dd4c-baea-4bb7-9a5d-9d05a637d44fn%40googlegroups.com.


FE_Nedelec.pdf
Description: Adobe PDF document


[deal.II] boundaries in dealii

2021-04-18 Thread Nikki Holtzer
I've recently implemented a 1d standard wave equation and have been messing 
around with different boundary conditions. For instance, I have run my code 
with 1 dirichlet and 1 neumann, 2 neumann, and 1 absorbing condition and 1 
neumann (as in step-24). When I run step-23 from the examples, which 
utilizes Dirichlet conditions, you can easily see the expected behavior of 
a phase flip after interaction from the boundary, i.e. as the source 
approaches the boundary it looks like a cap and after the boundary 
interaction it looks like a cup. However, in all 3 cases that I have tried, 
listed above, I am not observing the phase flip. No matter what boundary I 
change the left one to (dirichlet, neumann, absorbing) it seems to be 
giving me the same neumann behavior. Do you have any idea why this might be?

I have attached my code. If you were to run the code currently, it would 
run with a dirichlet condition on the left and a neumann on the right. In 
order to run with the absorbing boundary matrix , I comment out lines 
316-320 and line 611 as well as uncommenting the very end of lines 616 and 
645 which adds the boundary matrix to the computation. 

Thank you!!

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/eb188eba-ba35-4324-ac20-7bee8db4d810n%40googlegroups.com.
/* -
 *
 * This file aims to solve the standard wave equation u_tt-u_rr=0 and variants
 *
 * - */

#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 "helper.h"
#include 
#include 
using namespace dealii;


/* -- */
/* -- */


class Coefficients : public ParameterAcceptor
{
public:
  typedef std::complex value_type;
  static constexpr value_type imag{0., 1.};

  Coefficients()
  : ParameterAcceptor("A - Coefficients")
  {
R_0 = 0.0;
add_parameter("R_0", R_0, "inner radius of the computational domain");

R_1 = 10.0;
add_parameter("R_1", R_1, "outer radius of the computational domain");

A_param = 1.0;
add_parameter("A_param", A_param, "height of initial gaussian");

b_param = 5.0;
add_parameter("b_param", b_param, "position of peak of initial gaussian");

c_param = 2.0;
add_parameter("c_param", c_param, "wave speed of initial gaussian");

d_param = 1.0;
add_parameter("d_param", d_param, "standard deviation of initial gaussian");

f_param = 10.0;
add_parameter("f_param", f_param, "multiplicative constant of initial gaussian");

/* IC = A exp[-f{((x-b)+cT)^2/d^2}], initial_values1 = IC(0) */
   
initial_values1 = [&](double r) {
   return A_param*std::exp(-f_param * (r-b_param)*(r-b_param)/(pow(d_param,2.)));
};

//IC_T = d/dT(IC)(0)

IC_T = [&](double r) {
   return -(2.0 * c_param * f_param) / pow(d_param, 2.) * (r-b_param) * initial_values1(r);
};

boundary_values = [&](double r, double t) {
   if (r == R_0)
	 return 0.; 
	 //return t;
   if (r == R_1)
	 return 0.;
};

  }

  /* Publicly readable parameters: */
  double R_0;
  double R_1;
  double A_param;
  double b_param;
  double c_param;
  double d_param;
  double f_param;

  /* Publicly readable functions: */
  std::function initial_values1;
  std::function IC_T;
  std::function boundary_values;

};


/* -- */
/* -- */


template 
class Discretization : public ParameterAcceptor
{
public:
  static_assert(dim == 1, "Only implemented for 1D");

  Discretization(const Coefficients )
  : ParameterAcceptor("B - Discretization")
  , p_coefficients()
  {
ParameterAcceptor::parse_parameters_call_back.connect(
std::bind(::parse_parameters_callback, this));

refinement = 5;
add_parameter(
"refinement", refinement, "refinement of the spatial geometry");

order_mapping = 1;
add_parameter("order mapping", order_mapping, "order of the mapping");

order_finite_element = 1;
add_parameter("order finite element",
  order_finite_element,
  

[deal.II] Re: deal.II-based FSI matrix-free solver / looking for a post-doc

2021-04-18 Thread Peter Munch
Good luck at the defence! Hope it goes well!

Peter

On Sunday, 18 April 2021 at 00:44:44 UTC+2 mtwich...@gmail.com wrote:

> Dear all,
>  I have developed a parallel matrix-free monolithic FSI solver in the ALE 
> frame of reference. Exemplary results [1 , 2 
> ], more details in my PhD thesis [3] 
> 
> .
>
> *Highlights*:
>
> -> The proposed monolithic predictor-corrector time integration scheme 
> requires solving the generalized Stokes problem with discontinuous variable 
> coefficients. This has to be done one or two times per time step, depending 
> on the variant of the scheme.
>
> -> I have introduced a new multilevel preconditioner, for a 3D case the 
> time required to solve a system consisting of 421 million unknowns using 96 
> cores is 178 seconds.
>
> -> As demonstrated by the movie[1], I am able to solve 3D problems at 
> moderate Reynolds numbers (time step size of 0.01 did not result in any 
> stability issues with Re=2k). 
>
> ->The preconditioner from the thesis can still be greatly improved. I have 
> a proof of concept demonstrating the basic properties of a new idea of 
> smoother. The results look very promising and I'm curious how it would 
> work. Moreover, the algorithm is suitable for GPU that could result in 
> further speedup.
>  
> [1] 3D "Turek benchmark" at Re = 2000 (30M DoF)
> https://youtu.be/7TaC1dTfAQY
> [2] Water hammer in elastic tube (12M DoF)
> https://youtu.be/TRnlt_eC32Q
>
> Those results were obtained as a part of my PhD thesis. I have finally 
> managed to submit it in January, the defence will be held in 3 weeks 
> (expect reviews next week).
> I am looking for a post-doc (preferably deal.II-related). I'll be happy to 
> get an opportunity to continue my work.
>
> Best,
> Michał Wichrowski
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/c9fd51a4-cc3d-4fac-9da8-4f9a3fbaea66n%40googlegroups.com.