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

[Wolfgang Bangerth]

On 3/13/19 8:21 PM, Phạm Ngọc Kiên wrote:
> I am testing my codes for output the solution at a point in my numerical 
> model.
> I saw in the library 2 ways to do this task. The first one is to use 
> VectorTools::point_value() function and the second one is 
> fe_values.get_function_values().

These functions are really used in very different contexts. The first of these 
is when you want to evaluate a solution vector at some arbitrary point in the 
domain. As a consequence, you have to specify this point when you call the 
function.

The second function is used when you want to evaluate the function at very 
specific points, namely the quadrature points of the current cell. You can of 
course only use this if the point at which you want to evaluate the solution 
happens to be a quadrature point of the current cell, and not just any point 
anywhere. In return the second function is then vastly more efficient than the 
first.

Does this help?

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

[Phạm Ngọc Kiên]

Hi colleagues,
I am testing my codes for output the solution at a point in my numerical 
model.
I saw in the library 2 ways to do this task. The first one is to use 
VectorTools::point_value() function and the second one is 
fe_values.get_function_values().
For the latter function, I am using the code like this:

std::vector> temp_solution_re(1); 
std::vector> temp_solution_im(1);  //as my solution have both 
real and imaginary parts
fe_values[vector_re].get_function_values(solution, temp_solution_re);
fe_values[vector_im].get_function_values(solution, temp_solution_im);
for (unsigned int j = 0; j < dim; ++j) {
solution(j) = temp_solution_re[0][j];
solution(j+dim) = temp_solution_im[0][j];
}
Which returns the resutls are a vector of 3 component in x,y,z of my 
solution,and is the same as using
Vector values (dim+dim);
VectorTools::point_value(dof_handler,solution,point,values);

However, I checked that each tensor temp_solution_re[0], temp_solution_re[1], 
temp_solution_re[3] has 3 component in x,y,z direction.
Could you please tell me the meaning of these tensors when we get the solution 
by using fe_values.get_function_values() ?
Thank you very much.

Best regards,
Pham Ngoc Kien

-- 
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] ILU preconditioner and mesh-scale dependency?

[Jaekwang Kim]

Hi 

I had a trouble in solving scalar transport equation using DG elements. 

For example, I was solving for \lambda in the following equation where *u *is 
also a numerical solution

[image: Screen Shot 2019-03-13 at 9.08.17 PM.png]


The thing is I do not get any problem when my mesh file is scale of O(1) 
and I have tested my code works fine using the Method of Manufactured 
solution. 

Yet, when I tried to solve the governing equation at small scale 
mesh...(I.e. all (x,y) point in the original mesh-file was reduced by 
(0.01x,0.01y) 

I receive error like this 

An error occurred in line <129> of file 

 
in function

void dealii::SparseILU::initialize(const 
SparseMatrix &, const dealii::SparseILU::AdditionalData &) 
[number = double, somenumber = double]

The violated condition was: 

luval[ia[k]] != 0

Additional information: 

While computing the ILU decomposition, the algorithm found a zero pivot 
on the diagonal of row 34. This must stop the ILU algorithm because it 
means that the matrix for which you try to compute a decomposition is 
singular.


The problem occurs when the red line of the code is executed...

 template 

void StokesProblem::solve_transport ()

{

std::cout << "   -solve transport begins"<< std::endl;

SolverControl   solver_control (std::pow(10,6), 
system_rhs.block(2).l2_norm() * pow(10,-4));



unsigned int restart = 500;

SolverGMRES< Vector >::AdditionalData 
gmres_additional_data(restart+2);

SolverGMRES< Vector > solver(solver_control, 
gmres_additional_data);

  

//make preconditioner

SparseILU::AdditionalData additional_data(0,500); // (0 , 
additional diagonal terms)

std::cout << "B" << std::endl;

SparseILU preconditioner;



*preconditioner.initialize (system_matrix.block(2,2), 
additional_data);*



std::cout << "A" << std::endl;

solver.solve (system_matrix.block(2,2), solution.block(2), 
system_rhs.block(2), preconditioner);

// constraints.distribute (solution);

}




I wonder now ...

1. why smaller scale mesh only results such error,

2. If possible, is there any choice of better preconditionner for my system 
other than SparseILU to go around this problem? 


Thanks,

Regards, 

Jaekwang Kim 

-- 
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] Meeting some questions using constraints to apply boundary values

[Wolfgang Bangerth]

> I am coding for my graduation design, and hoping to achieve a coupling 
> geomechanics and fluid model.
> I use FE_RT for velocities and FE_Q for displacement of the rock.
> 
> In order to valid my model and programming, I use 1 D Terzaghi.
> But my boundary conditions appears have some problems.

I assume you then expect that the vertical velocity should be a 1d function 
independent of x and only dependent on y?

The picture you attach clearly does not satisfy this. I don't know why this 
isn't so, but here are a few things you can check:
* The vertical velocity you show is zero at the left and right sides of the 
domain. You say that you are applying "no flux" at boundaries 1, 2, 3 (which 
sides of the rectangle are these?) but that does not seem to be the case. I 
think that the function you are using, 
'VectorTools::project_boundary_values_div_conforming', should do this 
correctly, but you may want to look at a small testcase with just the 1x10 
mesh and output the constraints you get.

There is a function in DoFTools that also lets you output the physical 
location of every degree of freedom so that you can correlate what you get 
from the constraints object to where these DoFs are located.

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: About convergence using DG

[Wolfgang Bangerth]

On 3/13/19 9:18 AM, Apurva Tiwari wrote:
> 
> So far I'd been using an RK3 integrator that I had written myself.  To 
> rectify 
> the problem I faced (mentioned in the mail earlier),  I thought of, was using 
> the built-in explicit RK4 integrator TimeStepping::ExplicitRungeKutta 
> .
> Now, I have a function called "assemble_rhs" that uses MeshWorker to evaluate 
> the rhs of the explicit time dependent equation. Its signature is: 
> Step12::assemble_rhs (RHSIntegrator& rhs_integrator)
> (I'm modifying part of the step-12 tutorial program). But the evolve_one_step 
> function is of format:
> TimeStepping::ExplicitRungeKutta 
> <
>  
> VectorType >::evolve_one_time_step(   const std::function< 
> VectorType(const 
> double, const VectorType &)> &f,
>   double  t,
>   double  delta_t,
>   VectorType &y
>   )   
> 
> 
> To this, I tried to pass assemble_rhs as the first arg and ofcourse it failed 
> to compile. I've not been able to figure out how do I use my rhs function to 
> get to evaluate the next time-step. I'm probably making quite a fundamental 
> mistake, can you please guide me?

Yes, the function you have is not the function that works in the context of 
the time integrator.

You will want to take a look at the tutorial program that shows how to use the 
time integrator to see what kind of function can be passed as first argument.

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] How to apply first kind of boundary conditions to one component of a vector

[Wolfgang Bangerth]

On 3/10/19 7:39 AM, chenxi.edu1...@gmail.com wrote:
> 
> In order to check if this code is right, I want to use a model with 
> analytical 
> solution to verify it. And I need to make the displacement in y direction to 
> be zero.
> At first I use fe.component_mask(displacement),
> But this works for both displacements in x and y directions.
> 
> Then I try to construct a ComponentMask vector:
> I make the 2*dim+1 term to be false(which represent the constraints of y 
> direction)
> 
> However, I got this error:
> 
>   void
> dealii::VectorTools::{anonymous}::do_interpolate_boundary_values(const
> M_or_MC&, const DoFHandlerType&, const
> std::map*>&,
> std::map&, const dealii::ComponentMask&) [with
> int dim = 2; int spacedim = 2; number = double; DoFHandlerType =
> dealii::DoFHandler; M_or_MC = dealii::Mapping]
> 
> The violated condition was: 
> 
>      cell->get_fe().is_primitive (i)
> 
> Additional information: 
> 
>      This function can only deal with requested boundary values that
> correspond to primitive (scalar) base elements
> 
> 
> 
> Is there a good way to solve this problem?
> Thank you in advance!

You are trying to interpolate boundary values for components described by the 
Raviart-Thomas element. This element is not primitive -- see also
https://dealii.org/developer/doxygen/deal.II/DEALGlossary.html#GlossPrimitive
The problem with non-primitive elements is that one can not easily interpolate 
boundary values for them. But there are a number of other boundary value 
interpolation functions in namespace VectorTools that can be used for these 
kinds of elements.

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] Re: About convergence using DG

[Apurva Tiwari]

Hello!

So far I'd been using an RK3 integrator that I had written myself.  To
rectify the problem I faced (mentioned in the mail earlier),  I thought of,
was using the built-in explicit RK4 integrator
TimeStepping::ExplicitRungeKutta

.
Now, I have a function called "assemble_rhs" that uses MeshWorker to
evaluate the rhs of the explicit time dependent equation. Its signature is:
Step12::assemble_rhs (RHSIntegrator& rhs_integrator)
(I'm modifying part of the step-12 tutorial program). But the
evolve_one_step function is of format:
 TimeStepping::ExplicitRungeKutta
<
VectorType >::evolve_one_time_step ( const std::function< VectorType(const
double, const VectorType &)> &  f,
double  t,
double  delta_t,
VectorType &  y
)
To this, I tried to pass assemble_rhs as the first arg and ofcourse it
failed to compile. I've not been able to figure out how do I use my rhs
function to get to evaluate the next time-step. I'm probably making quite a
fundamental mistake, can you please guide me?

Thank You,
Apurva

On Tue, 12 Mar 2019 at 17:04, Apurva Tiwari  wrote:

> Hello!
> I'm trying to get appropriate convergence rates for a 2D scalar advection
> problem. The strategy is quite standard. I start with a coarse grid,
> successively refine it globally and generate the convergence table as done
> in step-7. As can be seen in the attached convergence tables, the values
> are far off. What all could have gone wrong and where do I start debugging
> my code from? What are the things to check?
> I'm using Q1, Q2 and Q3 elements with RK3 time-stepping.  The initial
> condition is an off-centered Gaussian which is subjected to a rotational
> velocity field.
>
> PFA: "error-global-order1.pdf "
>  "error-global-order2.pdf "
>  "error-global-order3.pdf "
>
> Thank You,
> Apurva Tiwari
>


-- 
Apurva Tiwari

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

[luca.heltai]

I believe the first point is a bug of the read_ucd function. 

I opened an issue for this (https://github.com/dealii/dealii/issues/7802)

The second point, however, is a feature. What should you do when two material 
ids meet? I think it is the responsability of the user to make sure that the 
correct manifold id is associated to the correct face.

Currently our read/write vtk functions (in the GridIn class) are the only ones 
allowing full treatment  material ids, boundary ids, and manifold ids.

I’m not sure UCD format allows you to specify both material and manifold ids 
for cells, and boundary and manifold is for faces. Maybe Andrea knows better?

L.

> On 8 Mar 2019, at 15:15, Nicola Giuliani  wrote:
> 
> PR 7775 has revealed some critical issues in the handling of manifold_id in 
> codim-1 domain.
> 
> -) 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.
> 
> -) 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?
> 

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