Re: [deal.II] Re: How to get the coordinates of a given degree of freedom on an edge?

2019-02-21 Thread Phạm Ngọc Kiên 
Hi,
The attachment is my small program that I created.
My fe system have 24 dofs for both real and imaginary of 12 edges of a cell.
No matter how I change the parameter "order" the number of dofs is
unchanged.
I think I have problem when constructing the object.
I would like to thank you very much.
 Best,
Kien

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

> On 2/20/19 7:01 PM, Phạm Ngọc Kiên wrote:
> > I tested in my codes using
> >
> > template
> > CSEM::CSEM()
> >  ://mapping(1),
> > dof_handler(triangulation),
> >  fe(FE_NedelecSZ(order),1,//(order), multiplicity
> > FE_NedelecSZ(order),1) {}
> >
> > When I changed the parameter order, I saw only the number of dofs was
> only of lowest order dofs no matter how big the order was
> >
> > I think that there should be a way to get higher order dofs, but I
> cannot find it now.
>
> Can you create a small program that demonstrates this issue?
>
> 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.
/* -
 *
 * Copyright (C) 2013 - 2018 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE.md at
 * the top level directory of deal.II.
 *
 * -


 */

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

#include 
#include 
#include 

#include 
#include 
#include 
#include 

#include 
#include 
#include 


#include 
#include 
#include 
#include 
#include 

#include 
#include 
#include 
#include 


#include 
#include 

#include 


#include 
#include 
#include 

#include 
#include 
#include 

using namespace dealii;

template
class EM {
public:
EM();


~EM();

void run();


private:
void grid_generation();

void setup_system();

void assemble_system();

void solve();



Triangulation triangulation;

DoFHandler dof_handler;
FESystem fe;

AffineConstraints constraints;

SparsityPattern sparsity_pattern;
SparseMatrix system_matrix;

Vector solution;
Vector system_rhs;

const unsigned int order = 3; //order of edge-based element

};

template
EM::EM()
:dof_handler(triangulation),
 fe(FE_NedelecSZ(order), 1,   //(order), multiplicity
FE_NedelecSZ(order), 1) {}

template
EM::~EM() {
dof_handler.clear();
}

template
void EM::grid_generation() {
GridGenerator::hyper_cube(triangulation, -4, 4);
}

template
void EM::setup_system() {
dof_handler.distribute_dofs(fe);

solution.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());

constraints.clear();
DoFTools::make_hanging_node_constraints(dof_handler, constraints);

VectorTools::project_boundary_values_curl_conforming_l2(dof_handler,
0,
ConstantFunction(0., fe.n_components()),
1,  
constraints);
VectorTools::project_boundary_values_curl_conforming_l2(dof_handler,
dim,
ConstantFunction(0., fe.n_components()),
1, 
constraints);
constraints.close();

[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-02-21 Thread chucui1016 
Dear Prof.Arndt,

Thank you very much! I understand what you say. It helps me a lot!

Best,
Chucui

在 2019年2月21日星期四 UTC+8下午7:31:20，Daniel Arndt写道：
>
> Chucui,
>
>
> For Question 1, I write a code to compute the L2 norm of (solution_1 - 
>> solution_2):
>> [...]
>>
>> Is that right?
>>
> You can just take the difference of the two variables and use 
> VectorTools::integrate_difference with a Functions::ZeroFunction as "exact 
> solution".
>  
>  
>
>>
>> For question 2, as the 2 solution vectors have different sizes, we don't 
>> have the same cells, so the cell loop cannot compute altogether:
>> typename DoFHandler::active_cell_iterator
>> cell_1 = dof_handler_1.begin_active(),
>> endc_1 = dof_handler_1.end(); 
>> typename DoFHandler::active_cell_iterator
>> cell_2 = dof_handler_2.begin_active();  
>> for (; cell_1!=endc_2; ++cell_1, ++cell_2)
>>   {}
>>
>> I wander how to solve this problem?
>>
> You need to interpolate the solution computed on the coarser mesh to the 
> finer mesh using SolutionTransfer(
> https://www.dealii.org/current/doxygen/deal.II/classSolutionTransfer.html, 
> explained e.g. in 
> https://www.dealii.org/current/doxygen/deal.II/step_26.html). Then, you 
> can take the difference just as in your first question.
>
> 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] deal.II at SIAM CSE (Spokane, WA)

2019-02-21 Thread Matthias Maier 
Dear all,

I had a quick chat with Wolfgang today and we were wondering who else is
attending the SIAM CSE conference in Spokane next week?

Also (in very German tradition) we decided that we absolutely have to
get an informal get-together of users and developers going!

Therefore, let's meet on Tuesday evening at 7pm at the "Onion Bar &
Grill" [1].

Please respond (either to the list or via private message) so that I can
get an estimate how many we will be.

Best,
Matthias

[1] No host - have a look at http://theonion.biz/ for menu and pricing

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

2019-02-21 Thread Rene Gassmoeller 
Hello everyone!

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

#7745: Fix FEFaceEvaluation with faces in non-standard orientation (proposed by 
kronbichler; merged) https://github.com/dealii/dealii/pull/7745

#7744: add type traits to be used internally with FEEValuation (proposed by 
davydden) https://github.com/dealii/dealii/pull/7744

#7743: Using rtrees in GridTools::compute_point_locations_try_all for cell 
search (proposed by GivAlz) https://github.com/dealii/dealii/pull/7743

#7742: Avoid ambiguous function declaration/variable initialization (proposed 
by masterleinad) https://github.com/dealii/dealii/pull/7742

#7741: Fix some doxygen errors (proposed by masterleinad; merged) 
https://github.com/dealii/dealii/pull/7741

#7740: Edits to the intro and results sections of step-61. (proposed by 
bangerth; merged) https://github.com/dealii/dealii/pull/7740

#7739: add random renumbering on level (proposed by tcclevenger; merged) 
https://github.com/dealii/dealii/pull/7739

#7738: Step-63 (proposed by tcclevenger) 
https://github.com/dealii/dealii/pull/7738

#7736: More HDF5 checks (proposed by masterleinad) 
https://github.com/dealii/dealii/pull/7736

#7735: Provide step-61 with the usual exception catching harness (proposed by 
masterleinad; merged) https://github.com/dealii/dealii/pull/7735

#7734: Remove redundant switch in MatrixFree::TaskInfo (proposed by 
masterleinad; merged) https://github.com/dealii/dealii/pull/7734

#7733: Correct MeshWorker::DoFInfo indices docstring (proposed by loganharbour; 
merged) https://github.com/dealii/dealii/pull/7733

#7732: allow zero padding in plate_with_a_hole (proposed by davydden; merged) 
https://github.com/dealii/dealii/pull/7732

#7731: [CI]: mark main status in github as failed (proposed by tjhei; merged) 
https://github.com/dealii/dealii/pull/7731

#7727: lac/linear_operator.h: Add rvalue-reference variant of inverse_operator 
(proposed by tamiko; merged) https://github.com/dealii/dealii/pull/7727

#7724: Adding cell_hint to GridTools::compute_point_locations documentation 
(proposed by GivAlz; merged) https://github.com/dealii/dealii/pull/7724

#7703: Revise iteration in PreconditionChebyshev (proposed by kronbichler; 
merged) https://github.com/dealii/dealii/pull/7703

#7578: Compute Point Location new function (proposed by GivAlz; merged) 
https://github.com/dealii/dealii/pull/7578

#6980: AD Helpers: Add helper for scalar functions (QP-level) (proposed by 
jppelteret; merged) https://github.com/dealii/dealii/pull/6980


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

#7737: doxygen doesn't run any more (opened) 
https://github.com/dealii/dealii/issues/7737


A list of all major changes since the last release can be found at 
https://www.dealii.org/developer/doxygen/deal.II/changes_after_8_5_0.html.


Thanks for being part of the community!


Let us know about questions, problems, bugs or just share your experience by 
writing to dealii@googlegroups.com, or by opening issues or pull requests at 
https://www.github.com/dealii/dealii.
Additional information can be found at https://www.dealii.org/.

-- 
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 get the coordinates of a given degree of freedom on an edge?

2019-02-21 Thread Wolfgang Bangerth 
On 2/20/19 7:01 PM, Phạm Ngọc Kiên wrote:
> I tested in my codes using
> 
> template
> CSEM::CSEM()
>  ://mapping(1),
> dof_handler(triangulation),
>  fe(FE_NedelecSZ(order),1,//(order), multiplicity
> FE_NedelecSZ(order),1) {}
> 
> When I changed the parameter order, I saw only the number of dofs was only of 
> lowest order dofs no matter how big the order was
> 
> I think that there should be a way to get higher order dofs, but I cannot 
> find it now.

Can you create a small program that demonstrates this issue?

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-02-21 Thread Wolfgang Bangerth 


> I am trying to solve a system of equations that do this:
> Stokes to solve for v_r and p_r for one fluid (viscous rock), I use these 
> solutions on the RHS of a Darcy type equation solved like step-20 for the 
> pressure p_f in the fluid in the domain. Using the 3 solutions, I update 
> another component (porosity - using a simple time dependent equation like 
> dphi/dt = rhs(pr-pf)) and put it back in the Stokes equation and so on
> 
> One of my conditions for the Darcy eqn is that for the top boundary, the 
> pressure there equals pr found previously.
> 
> Looking at the output, my values for p_f at the top boundary does not equal 
> the p_r just found. This is causing a massive issue as my porosity equation 
> needs pr-pf to be 0 (or at least very small) at the top boundary, and this 
> seems to be causing a real mess. I'm guessing because it's weakly imposed, 
> there might be issues.
> 
> Does anyone have any suggestions on what to do for this?
> 
> My base profile for the equations separately all work and have been verified 
> with MMS. With test values of the porosity, I have also verified that the 
> profiles separately give me what I expect (which then gives me pr-pf 
> monotonic 
> at least, but with what is happening at the top, this bit is messed up).

Jane -- the problem as you state it is too difficult to debug. Here's ho to 
approach this:

* Pick a problem for which you know the exact solution.
* In the Stokes equation, in all of the places where you currently use the 
Darcy solution, use the exact solution instead
* In the Darcy equation, in all of the places where you currently use the 
Stokes solution, use the exact solution instead

With this scheme, your iteration needs to converge in one step. If it doesn't, 
then you know where to look. If it does, revert the steps above one by one 
until the problem stops working.

I hope this helps!
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] Parallel Nonlinear Poisson code in DEALII

2019-02-21 Thread m boron 
Dear Jean-Paul,
Thank you very much for the quick response. Yes, it gives an error during
execution. I will also look at the Step-40 tutorial also.

Yours sincerely
Boron

On Thu, Feb 21, 2019 at 6:03 PM Jean-Paul Pelteret 
wrote:

> Dear Boron,
>
> As a new user, I’d like to welcome to the forum and deal.II!
>
> Am I correct that you only see this issue when executing your program in
> parallel? I think that the problem results from the way that you initialise
> the solution vector. That is, this line:
>
> present_solution.reinit(locally_owned_dofs, mpi_communicator);
>
>
> This is perfectly fine in many situations. However, in order to compute
> the solution values and gradients at quadrature points in a cell all
> degrees of freedom belonging to that cell must be accessible from the
> solution vector. Since this vector “ present_solution.reinit" does not
> contain the ghost DoFs, there will eventually be a cell that has its DoFs
> split between processes and not all values will be immediately accessible
> on a process. The way to fix this issue would be to initialise an
> intermediate vector that also holds ghost values, copy the distributed
> solution to the ghosted one, and then proceed with all cell-based
> calculations using the ghosted vector. That would be done something like
> this:
>
>   // At the beginning of assemble_system()
>   PETScWrappers::MPI::Vector present_solution_ghosted;
>   present_solution_ghosted.reinit(locally_owned_dofs,
> locally_relevant_dofs , mpi_communicator);
>   …
>   // On a cell
>fe_values.get_function_values(present_solution_ghosted, old_solution);
>
> A tutorial that shows how to get the locally_relevant_dofs is step-40
> .
>
> I hope that this helps you solve the problem.
> Best,
> Jean-Paul
>
> On 21 Feb 2019, at 11:02, mboron1...@gmail.com wrote:
>
> Dear all,
> I am Boron. I am a new to DEALII. I am currently trying to write a
> parallel code in DEALII for solving nonlinear Poisson's equation. The file
> is also attahed below. My doubt is "How do we pass history variable while
> constructing the cell_matrix?"
> A code snippet is (Line No 211-225) :
>
> for (; cell!=endc; ++cell)
> {
> if (cell->subdomain_id() == this_mpi_process)
> {
> cell_matrix = 0;
> cell_rhs = 0;
>
> fe_values.reinit (cell);
>
> fe_values.get_function_values(present_solution,
> old_solution);
> fe_values.get_function_gradients(present_solution,
> old_solution_gradients);
> for (unsigned int q_point=0; q_point {
> // BUILD ELEMENTAL CELL MATRIX @ EACH GAUSS POINT
> }
>
> In the above code snippet, the line
> 'fe_values.get_function_values(present_solution, old_solution);  ' throws
> an error. Is there a way to pass only the vectors relevant to the
> corresponding subdomains in DEALII?
>
>
> --
> 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.
> http://nonlinearpoissonparallel.cc>>
>
>
> --
> 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] Parallel Nonlinear Poisson code in DEALII

2019-02-21 Thread Jean-Paul Pelteret 
Dear Boron,

As a new user, I’d like to welcome to the forum and deal.II!

Am I correct that you only see this issue when executing your program in 
parallel? I think that the problem results from the way that you initialise the 
solution vector. That is, this line:

> present_solution.reinit(locally_owned_dofs, mpi_communicator);


This is perfectly fine in many situations. However, in order to compute the 
solution values and gradients at quadrature points in a cell all degrees of 
freedom belonging to that cell must be accessible from the solution vector. 
Since this vector “ present_solution.reinit" does not contain the ghost DoFs, 
there will eventually be a cell that has its DoFs split between processes and 
not all values will be immediately accessible on a process. The way to fix this 
issue would be to initialise an intermediate vector that also holds ghost 
values, copy the distributed solution to the ghosted one, and then proceed with 
all cell-based calculations using the ghosted vector. That would be done 
something like this:

  // At the beginning of assemble_system()
  PETScWrappers::MPI::Vector present_solution_ghosted;
  present_solution_ghosted.reinit(locally_owned_dofs, locally_relevant_dofs , 
mpi_communicator);
  …
  // On a cell
   fe_values.get_function_values(present_solution_ghosted, old_solution);  

A tutorial that shows how to get the locally_relevant_dofs is step-40 
.

I hope that this helps you solve the problem.
Best,
Jean-Paul  

> On 21 Feb 2019, at 11:02, mboron1...@gmail.com wrote:
> 
> Dear all,
> I am Boron. I am a new to DEALII. I am currently trying to write a parallel 
> code in DEALII for solving nonlinear Poisson's equation. The file is also 
> attahed below. My doubt is "How do we pass history variable while 
> constructing the cell_matrix?"
> A code snippet is (Line No 211-225) :
> 
> for (; cell!=endc; ++cell)
> {
> if (cell->subdomain_id() == this_mpi_process)
> {
> cell_matrix = 0; 
> cell_rhs = 0;
> 
> fe_values.reinit (cell);
> 
> fe_values.get_function_values(present_solution, old_solution);
> fe_values.get_function_gradients(present_solution, 
> old_solution_gradients);
> for (unsigned int q_point=0; q_point {
> // BUILD ELEMENTAL CELL MATRIX @ EACH GAUSS POINT
> }
> 
> In the above code snippet, the line 
> 'fe_values.get_function_values(present_solution, old_solution);  ' throws an 
> error. Is there a way to pass only the vectors relevant to the corresponding 
> subdomains in DEALII?
> 
> 
> -- 
> 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: 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-02-21 Thread Daniel Arndt 
Chucui,


For Question 1, I write a code to compute the L2 norm of (solution_1 - 
> solution_2):
> [...]
>
> Is that right?
>
You can just take the difference of the two variables and use 
VectorTools::integrate_difference with a Functions::ZeroFunction as "exact 
solution".
 
 

>
> For question 2, as the 2 solution vectors have different sizes, we don't 
> have the same cells, so the cell loop cannot compute altogether:
> typename DoFHandler::active_cell_iterator
> cell_1 = dof_handler_1.begin_active(),
> endc_1 = dof_handler_1.end(); 
> typename DoFHandler::active_cell_iterator
> cell_2 = dof_handler_2.begin_active();  
> for (; cell_1!=endc_2; ++cell_1, ++cell_2)
>   {}
>
> I wander how to solve this problem?
>
You need to interpolate the solution computed on the coarser mesh to the 
finer mesh using 
SolutionTransfer(https://www.dealii.org/current/doxygen/deal.II/classSolutionTransfer.html,
 
explained e.g. in 
https://www.dealii.org/current/doxygen/deal.II/step_26.html). Then, you can 
take the difference just as in your first question.

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

2019-02-21 Thread Daniel Arndt 
Jane, 

can you write down the discretized weak formulation you using?
Is there any paper showing well-posedness and suitable error estimates for 
that discretization?

Best,
Daniel

Am Mittwoch, 6. Februar 2019 17:22:28 UTC+1 schrieb jane...@jandj-ltd.com:
>
> Hi all,
>
> I am trying to solve a system of equations that do this:
> Stokes to solve for v_r and p_r for one fluid (viscous rock), I use these 
> solutions on the RHS of a Darcy type equation solved like step-20 for the 
> pressure p_f in the fluid in the domain. Using the 3 solutions, I update 
> another component (porosity - using a simple time dependent equation like 
> dphi/dt = rhs(pr-pf)) and put it back in the Stokes equation and so on 
>
> One of my conditions for the Darcy eqn is that for the top boundary, the 
> pressure there equals pr found previously. 
>
> Looking at the output, my values for p_f at the top boundary does not 
> equal the p_r just found. This is causing a massive issue as my porosity 
> equation needs pr-pf to be 0 (or at least very small) at the top boundary, 
> and this seems to be causing a real mess. I'm guessing because it's weakly 
> imposed, there might be issues. 
>
> Does anyone have any suggestions on what to do for this? 
>
> My base profile for the equations separately all work and have been 
> verified with MMS. With test values of the porosity, I have also verified 
> that the profiles separately give me what I expect (which then gives me 
> pr-pf monotonic at least, but with what is happening at the top, this bit 
> is messed up).
>
>
> For the boundary condition, I am doing:
>
> for (unsigned int face_no=0;
>
> face_no::faces_per_cell;
>
> ++face_no)
>
> if (cell->face(face_no)->at_boundary()
>
> &&
>
> (cell->face(face_no)->boundary_id() == 1)) // top has boundary id 1
>
> {
>
> fe_face_values_pf.reinit (cell, face_no);
>
> fe_face_values_rock.reinit (vr_cell, face_no);
>
>
> fe_face_values_rock[pressure].get_function_values (solution_rock, 
> pr_boundary_values);
>
>
> // DIRICHLET CONDITION FOR TOP. pf = pr at top boundary
>
> for (unsigned int q=0; q
> {
>
> for (unsigned int i=0; i
> {
>
> local_rhs(i) += -( fe_face_values_pf[velocities].value (i, q) *
>
> fe_face_values_pf.normal_vector(q) *
>
> pr_boundary_values[q] *
>
> fe_face_values_pf.JxW(q));
>
> }
>
> }
>
> }
>
>

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

2019-02-21 Thread jane . lee 
Hi all,
Does anyone have any suggestions on this? I'm still struggling to get the 
values to equal more precisely.

Thanks


On Wednesday, February 6, 2019 at 7:22:28 PM UTC+3, jane...@jandj-ltd.com 
wrote:
>
> Hi all,
>
> I am trying to solve a system of equations that do this:
> Stokes to solve for v_r and p_r for one fluid (viscous rock), I use these 
> solutions on the RHS of a Darcy type equation solved like step-20 for the 
> pressure p_f in the fluid in the domain. Using the 3 solutions, I update 
> another component (porosity - using a simple time dependent equation like 
> dphi/dt = rhs(pr-pf)) and put it back in the Stokes equation and so on 
>
> One of my conditions for the Darcy eqn is that for the top boundary, the 
> pressure there equals pr found previously. 
>
> Looking at the output, my values for p_f at the top boundary does not 
> equal the p_r just found. This is causing a massive issue as my porosity 
> equation needs pr-pf to be 0 (or at least very small) at the top boundary, 
> and this seems to be causing a real mess. I'm guessing because it's weakly 
> imposed, there might be issues. 
>
> Does anyone have any suggestions on what to do for this? 
>
> My base profile for the equations separately all work and have been 
> verified with MMS. With test values of the porosity, I have also verified 
> that the profiles separately give me what I expect (which then gives me 
> pr-pf monotonic at least, but with what is happening at the top, this bit 
> is messed up).
>
>
> For the boundary condition, I am doing:
>
> for (unsigned int face_no=0;
> face_no::faces_per_cell;
> ++face_no)
> if (cell->face(face_no)->at_boundary()
> &&
> (cell->face(face_no)->boundary_id() == 1)) // top has boundary id 1
> {
> fe_face_values_pf.reinit (cell, face_no);
> fe_face_values_rock.reinit (vr_cell, face_no);
> fe_face_values_rock[pressure].get_function_values (solution_rock, 
> pr_boundary_values);
> // DIRICHLET CONDITION FOR TOP. pf = pr at top boundary
> for (unsigned int q=0; q {
> for (unsigned int i=0; i {
> local_rhs(i) += -( fe_face_values_pf[velocities].value (i, q) *
> fe_face_values_pf.normal_vector(q) *
> pr_boundary_values[q] *
> fe_face_values_pf.JxW(q));
> }
> }
> }
>
>

-- 
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] Parallel Nonlinear Poisson code in DEALII

2019-02-21 Thread mboron1982 
Dear all,
I am Boron. I am a new to DEALII. I am currently trying to write a parallel 
code in DEALII for solving nonlinear Poisson's equation. The file is also 
attahed below. My doubt is "How do we pass history variable while 
constructing the cell_matrix?"
A code snippet is (Line No 211-225) :

for (; cell!=endc; ++cell)
{
if (cell->subdomain_id() == this_mpi_process)
{
cell_matrix = 0; 
cell_rhs = 0;

fe_values.reinit (cell);

fe_values.get_function_values(present_solution, 
old_solution);
fe_values.get_function_gradients(present_solution, 
old_solution_gradients);
for (unsigned int q_point=0; q_pointhttp://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.
/* PROGRAM TO SOLVE NON-LINEAR POISSON IN PARALLEL USING PETSC*/
// NABLA DOT ((1+U) (NABLA U)) = F
// U =sin(3x) * sin(3y) in (0,1)^2

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

#include 
#include 
#include 
#include 

#include 
#include 
#include 

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

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

#include 
#include 
#include 
#include 

#include 
#include 
#include  
#include 

#include 

namespace Non_Linear_Poisson
{
	using namespace std;
	using namespace dealii;

	template 
	class NLPoisson
	{
	public:
		NLPoisson();
		~NLPoisson();
		void run();

	private:
		void setup_system(const bool initial_step);
		void assemble_system();
		unsigned int solve();
		void refine_mesh();
		void set_boundary_values();
		double compute_residual(const double alpha) const;
		double determine_step_length() const;
	
		MPI_Comm mpi_communicator;
		
	const unsigned int n_mpi_processes;
	const unsigned int this_mpi_process;
	
	ConditionalOStream pcout;
		
		Triangulation 			triangulation;
		DoFHandler			dof_handler; 
	
		FE_Q	 fe;
	
		ConstraintMatrix	 		hanging_node_constraints;
	
	PETScWrappers::MPI::SparseMatrix 		system_matrix;
		PETScWrappers::MPI::Vectorpresent_solution;
		PETScWrappers::MPI::Vectornewton_update;
		PETScWrappers::MPI::Vectorsystem_rhs;
	};

	template 
	class BoundaryValues : public Function
	{
	public:
		BoundaryValues(): Function () {}
		
		virtual double value(const Point , const unsigned int component = 0) const;
	};
	
	template
	double BoundaryValues::value(const Point , const unsigned int /*component*/) const
	{
		return sin(3*p[0]) * sin(3*p[1]);
	}
	
	template 
	class RightHandSide : public Function
	{
	public:
		RightHandSide(): Function () {}
		
		virtual double value(const Point , const unsigned int component = 0) const;
	};
	
	template
	double RightHandSide::value(const Point , const unsigned int /*component*/) const
	{		
		return 18 * sin(3*p[0]) * sin(3*p[1])* (sin(3*p[0]) * sin(3*p[1]) + 1) 	- 
9*cos(3*p[1])*cos(3*p[1]) * sin(3*p[0])*sin(3*p[0]) 			- 
9*cos(3*p[0])*cos(3*p[0]) * sin(3*p[1])*sin(3*p[1]);
	}
	

	template
	NLPoisson::NLPoisson() 
		: 
		mpi_communicator (MPI_COMM_WORLD),
		n_mpi_processes (Utilities::MPI::n_mpi_processes(mpi_communicator)),
	this_mpi_process (Utilities::MPI::this_mpi_process(mpi_communicator)),
	pcout (std::cout, (this_mpi_process == 0)),
		dof_handler(triangulation), fe(1)
	{}

	template
	NLPoisson::~NLPoisson() 
	{
		dof_handler.clear();	
	}

	template 
	void NLPoisson::setup_system(const bool initial_step)
	{
		if (initial_step)
		{
			pcout<<"	inside setup 1"< locally_owned_dofs_per_proc = DoFTools::locally_owned_dofs_per_subdomain(dof_handler);
		const IndexSet locally_owned_dofs = locally_owned_dofs_per_proc[this_mpi_process];
			present_solution.reinit(locally_owned_dofs, mpi_communicator);

			
		}
pcout<<"	inside setup 2"< locally_owned_dofs_per_proc = DoFTools::locally_owned_dofs_per_subdomain(dof_handler);
	const IndexSet locally_owned_dofs = locally_owned_dofs_per_proc[this_mpi_process];
	
		system_matrix.reinit(locally_owned_dofs, locally_owned_dofs, dsp, mpi_communicator);
		newton_update.reinit (locally_owned_dofs, mpi_communicator);
	system_rhs.reinit (locally_owned_dofs, mpi_communicator);

	}

	template 
  	void NLPoisson::assemble_system ()
  	{
  		pcout<<"	inside assembly"<  quadrature_formula(3);
		
	FEValues fe_values (fe, quadrature_formula,
 update_values | update_gradients |
 update_quadrature_points | update_JxW_values);

	const unsigned int   dofs_per_cell = fe.dofs_per_cell;
	const unsigned

Re: [deal.II] 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-02-21 Thread chucui1016 
Dear Jean-Paul,

Thank you very much for your quick reply! I will read it now! 

Best,
Chucui

在 2019年2月21日星期四 UTC+8下午5:17:17，Jean-Paul Pelteret写道：
>
> Hi Chucui,
>
> I can only offer you a very brief reply right now: You might want to look 
> into the ConvergenceTable 
>  
> class 
> to help you with (1). The introduction to the class directs you to some 
> tutorials that indicate its use.
>
> Best,
> Jean-Paul
>
> On 21 Feb 2019, at 10:11, chucu...@gmail.com  wrote:
>
> Hi all,
>
> I have 2 questions of computing convergence rate of L2 norm of error:
>
> 1. If I don't have exact solution, does deal.ii have functions to compute 
> convergence rate directly? (Like VectorTools::integrate_difference, which 
> deals with the problem with exact solution  )  
>   
> 2. If  I want to compute convergence rate of error of solutions with 
> different globally refinement (u_{h}-u_{h/2}, u_{h/2}-u_{h/4}, and so on ), 
> does deal.ii have  functions to compute convergence rate directly?
> if not, how to deal with computing u_{h}-u_{h/2}, where the 2 numerical 
> solution have different size of vector?
>
> Thanks in advance!
>
> Best,
> Chucui
>
> -- 
> 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.

[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-02-21 Thread chucui1016 
Hi, all,

If there is no function to compute convergence above, I need to write the 
code by myself.

For Question 1, I write a code to compute the L2 norm of (solution_1 - 
solution_2):

template 
  double StokesProblem::norm_compute (const Vector solu_1,
  const Vector solu_2)
  {
QGauss   quadrature_formula(degree+2);

FEValues fe_values (fe, quadrature_formula,
 update_values|
 update_quadrature_points  |
 update_JxW_values |
 update_gradients);
  
const unsigned int   dofs_per_cell   = fe.dofs_per_cell;

const unsigned int   n_q_points  = quadrature_formula.size();
FullMatrix   local_con_rate_matrix (dofs_per_cell, 
dofs_per_cell);
std::vector local_dof_indices (dofs_per_cell); 

std::vector solu_1_values(n_q_points);
std::vector solu_2_values(n_q_points);

double norm_12 = 0, norm_24 = 1,norm_12_square = 0, norm_24_square = 1, 
rate_124 = 0, con_rate = 0; 
typename DoFHandler::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();   
for (; cell!=endc; ++cell)
  {
fe_values.reinit (celli);  
fe_values.get_function_values (solu_1, solu_1_values);
fe_values.get_function_values (solu_2, solu_2_values);
for (unsigned int q=0; q::active_cell_iterator
cell_1 = dof_handler_1.begin_active(),
endc_1 = dof_handler_1.end(); 
typename DoFHandler::active_cell_iterator
cell_2 = dof_handler_2.begin_active();  
for (; cell_1!=endc_2; ++cell_1, ++cell_2)
  {}

I wander how to solve this problem?

Thanks in advance!
Best,
Chucui

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

2019-02-21 Thread Jean-Paul Pelteret 
Hi Chucui,

I can only offer you a very brief reply right now: You might want to look into 
the ConvergenceTable 
 class 
to help you with (1). The introduction to the class directs you to some 
tutorials that indicate its use.

Best,
Jean-Paul

> On 21 Feb 2019, at 10:11, chucui1...@gmail.com wrote:
> 
> Hi all,
> 
> I have 2 questions of computing convergence rate of L2 norm of error:
> 
> 1. If I don't have exact solution, does deal.ii have functions to compute 
> convergence rate directly? (Like VectorTools::integrate_difference, which 
> deals with the problem with exact solution  )  
>   
> 2. If  I want to compute convergence rate of error of solutions with 
> different globally refinement (u_{h}-u_{h/2}, u_{h/2}-u_{h/4}, and so on ), 
> does deal.ii have  functions to compute convergence rate directly?
> if not, how to deal with computing u_{h}-u_{h/2}, where the 2 numerical 
> solution have different size of vector?
> 
> Thanks in advance!
> 
> Best,
> Chucui
> 
> -- 
> 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] 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-02-21 Thread chucui1016 
Hi all,

I have 2 questions of computing convergence rate of L2 norm of error:

1. If I don't have exact solution, does deal.ii have functions to compute 
convergence rate directly? (Like VectorTools::integrate_difference, which 
deals with the problem with exact solution  )  
  
2. If  I want to compute convergence rate of error of solutions with 
different globally refinement (u_{h}-u_{h/2}, u_{h/2}-u_{h/4}, and so on ), 
does deal.ii have  functions to compute convergence rate directly?
if not, how to deal with computing u_{h}-u_{h/2}, where the 2 numerical 
solution have different size of vector?

Thanks in advance!

Best,
Chucui

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