Re: [deal.II] Reuse of the system_matrix in the Picard iteration for Inhomogenous Dirichlet boundary conditions

2019-04-26 Thread Wolfgang Bangerth
On 4/26/19 8:04 PM, hanks0...@gmail.com wrote:
> 
> 
> The several components in the matrix A_ij change to 0 or 1, and these changed 
> values are used in assembling the system_rhs.
> (if it is wrong please let me know what is the reason)
> 
> So, my question is reassembling the system_matrix is necessary...?

Not really. You will want to use something like what is done in step-26 where 
we store the matrix (actually, matrices) without boundary values applied, and 
in each time step copy them into system_matrix and apply boundary values. I 
suspect that this is what you also want to do in your case.

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: Integrating a deal.II-specific function in a NOX-Interface for MPI-threads

2019-04-26 Thread Wolfgang Bangerth
On 4/26/19 2:28 PM, 'Maxi Miller' via deal.II User Group wrote:
> I think the current version is the smallest version I can get. It will work 
> for a single thread, but not for two or more threads.

I tried to run this on my machine, but I don't have NOX as part of my Trilinos 
installation. That might have to wait for next week.

But I think there is still room for making the program shorter. The point is 
that the parts of the program that run before the error happens do not 
actually have to make sense. For example, could you replace your own boundary 
values class by ZeroFunction? Instead of assembling something real, could you 
just use the identity matrix on each cell?

Other things you don't need: Timers, convergence tables, etc. Which of the 
parameters you set in solve_for_NOX() are necessary? (Necessary to reproduce 
the bug; I know they're necessary for your application, but that's irrelevant 
here.) Could you replace building the residual by just putting random stuff in 
the vector?

Similarly, the computeJacobian and following functions really just check 
conditions, but when you run the program to illustrate the bug you're chasing, 
do these checks trigger? If not, remove them, then remove the calls to these 
functions (because they don't do anything any more), then remove the whole 
function.

I think there's still room to make this program smaller by at least a factor 
of 2 or 3 or 4 :-)

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] Disordered mesh elements after several refinmnets

2019-04-26 Thread Wolfgang Bangerth
On 4/26/19 8:35 PM, ofcrashb...@gmail.com wrote:
> 
> I am solving Laplace equation near end of line where field is changing like 
> sqrt(r), thats why I dense mesh size near tip, and also use adaptive mesh 
> refinement to make it dense even more.
> 
> But resulted mesh is very disordered(even simulation still works).
> I can't understand where is the problem: in Paraview, in mesh 
> generator(Triangle+Tethex) or in Deal.ii? And what are others limitation of 
> mesh size, cos I need those to be very small near tip.

What do the pictures you attached show? It looks like some of the cells are 
triangles, not quadrilaterals. Is this a 3d mesh? A 2d mesh?

Best
  W.


-- 

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

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


Re: [deal.II] Imposing the Neumann BC using right_hand_side() function

2019-04-26 Thread Wolfgang Bangerth
On 4/26/19 3:42 AM, Muhammad Mashhood wrote:
> Thank you Prof. Bangerth for the correction that was quit helping. Currently 
> the code started running but unfortunately the function is returning with the 
> zero "tmp" vector. Any possible correction?

Hm, I wonder if it is even implemented in the 1d case. Looking into the 
deal.II source files, I can't see how your program actually runs without 
producing an error. Are you using debug mode?

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] Step 12, rewritten using ScratchData and CopyData

2019-04-26 Thread 'Maxi Miller' via deal.II User Group
I tried to rewrite step-12 using the recently introduced classes 
ScratchData and CopyData. In the first iteration, the results are similar, 
in the second iteration the new version gives wrong results, though. I 
compared the generated matrices and right hand side vectors, and they were 
identical, but after the first mesh refinement the rewritten version 
results in wrong solutions. Thus I was wondering if I forgot to initialize 
something, or if I forgot to loop over/add some cells after the refinement?
Thanks!

-- 
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) 2009 - 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.
 *
 * -

 *
 * Author: Guido Kanschat, Texas A University, 2009
 */


// The first few files have already been covered in previous examples and will
// thus not be further commented on:
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
// Here the discontinuous finite elements are defined. They are used in the same
// way as all other finite elements, though -- as you have seen in previous
// tutorial programs -- there isn't much user interaction with finite element
// classes at all: they are passed to DoFHandler and
// FEValues objects, and that is about it.
#include 
// We are going to use the simplest possible solver, called Richardson
// iteration, that represents a simple defect correction. This, in combination
// with a block SSOR preconditioner (defined in precondition_block.h), that
// uses the special block matrix structure of system matrices arising from DG
// discretizations.
#include 
#include 
// We are going to use gradients as refinement indicator.
#include 

// Here come the new include files for using the MeshWorker framework. The first
// contains the class MeshWorker::DoFInfo, which provides local integrators with
// a mapping between local and global degrees of freedom. It stores the results
// of local integrals as well in its base class MeshWorker::LocalResults.
// In the second of these files, we find an object of type
// MeshWorker::IntegrationInfo, which is mostly a wrapper around a group of
// FEValues objects. The file meshworker/simple.h contains classes
// assembling locally integrated data into a global system containing only a
// single matrix. Finally, we will need the file that runs the loop over all
// mesh cells and faces.
#include 
#include 
#include 
#include 

// Like in all programs, we finish this section by including the needed C++
// headers and declaring we want to use objects in the dealii namespace without
// prefix.
#include 
#include 


namespace Step12
{
  using namespace dealii;

  // @sect3{Equation data}
  //
  // First, we define a class describing the inhomogeneous boundary data. Since
  // only its values are used, we implement value_list(), but leave all other
  // functions of Function undefined.
  template 
  class BoundaryValues : public Function
  {
  public:
BoundaryValues() = default;
virtual void value_list(const std::vector> ,
std::vector &  values,
const unsigned int component = 0) const override;
  };

  // Given the flow direction, the inflow boundary of the unit square $[0,1]^2$
  // are the right and the lower boundaries. We prescribe discontinuous boundary
  // values 1 and 0 on the x-axis and value 0 on the right boundary. The values
  // of this function on the outflow boundaries will not be used within the DG
  // scheme.
  template 
  void BoundaryValues::value_list(const std::vector> ,
   std::vector &  values,
   const unsigned int component) const
  {
(void)component;
AssertIndexRange(component, 1);
Assert(values.size() == points.size(),
   ExcDimensionMismatch(values.size(), points.size()));

for (unsigned int i = 0; i < values.size(); ++i)
  {
if 

Re: [deal.II] Re: Integrating a deal.II-specific function in a NOX-Interface for MPI-threads

2019-04-26 Thread 'Maxi Miller' via deal.II User Group
I think the current version is the smallest version I can get. It will work 
for a single thread, but not for two or more threads.
One observation I made was that if I not only copy the locally owned 
values, but also the locally relevant values from the deal.II-vector to the 
Epetra_Vector, the solver still does not converge, but the result is 
correct. When changing line 159 from locally_owned_elements to 
locally_relevant_elements, I get an access error, though, so there must be 
a different way to achieve that.

Am Freitag, 26. April 2019 17:00:58 UTC+2 schrieb Wolfgang Bangerth:
>
> On 4/26/19 1:59 AM, 'Maxi Miller' via deal.II User Group wrote: 
> > I tried to reduce it as much as possible while still being able to run, 
> > the result is attached. It still has ~600 lines, but I am not sure if I 
> > should remove things like the output-function? Furthermore the NOX-logic 
> > takes quite a lot of lines, and I would like to keep the program as 
> > usable as possible while removing unneccessary lines. Is that example 
> > usable? 
>
> Remove everything you can. If you get an error at one point, no output 
> will be created -- so remove the output function. If the error happens 
> before the linear system is solved, remove the code that computes the 
> entries of the matrix. Then remove the matrix object itself. If you are 
> currently solving a linear system before the error happens, try what 
> happens if you just comment out the solution step -- and if the error 
> still happens (because it probably doesn't depend on the actual values 
> of the vector), then remove the solution step and the assembly of the 
> matrix and the matrix. If you need a DoFHandler to set up the vector, 
> output the index sets and sizes you use for the processors involved and 
> build these index sets by hand instead -- then remove the DoFHandler and 
> the finite element and everything else related to this. Continue this 
> style of work until you really do have a small testcase :-) 
>
> Cheers 
>   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.
//Nox include files
#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 

// NOX Objects
#include "NOX.H"
#include "NOX_Epetra.H"

// Trilinos Objects
#ifdef HAVE_MPI
#include "Epetra_MpiComm.h"
#else
#include "Epetra_SerialComm.h"
#endif
#include "Epetra_Vector.h"
#include "Epetra_Operator.h"
#include "Epetra_RowMatrix.h"
#include "NOX_Epetra_Interface_Required.H" // base class
#include "NOX_Epetra_Interface_Jacobian.H" // base class
#include "NOX_Epetra_Interface_Preconditioner.H" // base class
#include "Epetra_CrsMatrix.h"
#include "Epetra_Map.h"
#include "Epetra_LinearProblem.h"
#include "AztecOO.h"
#include "Teuchos_StandardCatchMacros.hpp"

// User's application specific files
//#include "problem_interface.h" // Interface file to NOX
//#include "finite_element_problem.h"

// Required for reading and writing parameter lists from xml format
// Configure Trilinos with --enable-teuchos-extended
#ifdef HAVE_TEUCHOS_EXTENDED
#include "Teuchos_XMLParameterListHelpers.hpp"
#include 
#endif

// This is C++ ...
#include 
#include 

using namespace dealii;

template 
class BoundaryValues : public dealii::Function
{
public:
BoundaryValues(const size_t n_components)
: dealii::Function(1), n_components(n_components)
{}

virtual double value (const dealii::Point   ,
  const unsigned int  component = 0) const;

virtual void vector_value(const dealii::Point , dealii::Vector ) const;
private:
const size_t n_components;
};


template 
double BoundaryValues::value (const dealii::Point ,
   const unsigned int) const
{
const double x = p[0];
const double y = p[1];
return sin(M_PI * x) * sin(M_PI * y);
}

template 
void BoundaryValues::vector_value(const dealii::Point , dealii::Vector ) const
{
for(size_t i = 0; i < value.size(); ++i)
value[i] = BoundaryValues::value(p, i);
}

template 
class ProblemInterface : public 

Re: [deal.II] Diffusion equation with DG fails with distorted mesh

2019-04-26 Thread Wolfgang Bangerth
On 4/26/19 9:48 AM, Gary Uppal wrote:
> 
> Thanks for the reply. I get problems on less distorted meshes as well. 
> For example, in the mesh with a hole shown below.  I eventually need 
> multiple holes and an advection term. Actually, the advection seems ok 
> and I mainly get issues when the diffusion is large. If the distorted 
> mesh is the issue, is there a better way to construct the geometry?

No, your mesh is reasonable. The question is not whether the error is 
large on cells that are not squares (it may be -- the finite element 
theory doesn't ever guarantee that the error is *small* on any given 
cell of any given mesh), but whether you get the *correct convergence 
order*. So looking at the solution on a single mesh is not a useful 
strategy -- the question is how the error decays as you refine the mesh.

Best
  W.


-- 

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

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


Re: [deal.II] Diffusion equation with DG fails with distorted mesh

2019-04-26 Thread Wolfgang Bangerth
On 4/26/19 2:50 AM, Doug wrote:
> 
> I recover optimal convergence orders (p+1) for straight meshes when 
> refining in h. I recover optimal orders (p+1) for a sine-transformed 
> grid with some skewed angles up to around 45 degrees (as seen attached).
> 
> However, I lose my optimal (p+1) order as soon as I apply any small 
> amount of random distortion to my final refined grid (see attached). 
> Note that the distortion factor of 0.15 is used, but this loss of order 
> also occurs for small distortion factors down to 0.01. Instead of (p+1), 
> I may recover (p) or even (p-1). This only applies when the final grid 
> is randomly distorted. If I start from a coarse-ish distorted grid (say 
> 100 cells), and refine globally such that its children aren't distorted, 
> I then recover my optimal orders again.

I don't know enough about DG theory to tell what this is from. But this 
is true:
* for your example with the sine-transformed domain, if you refine
   the mesh sufficiently many times, each cell will be getting closer and
   closer to a parallelogram
* the same is the case for your initially-distorted then refined mesh.
As a matter of fact, that's a general theorem: Start with some mesh and 
refine it sufficiently often, and all cells with tend to parallelograms.

Why does this matter? Because the mapping from the reference cell to a 
parallelogram is linear, and consequently the derivative of the mapping 
(which shows up in the convergence proofs of the finite element method) 
will be a constant on every cell in the limit of h->0.

On the other hand, if you refine a mesh *and then distort them 
randomly*, then the mapping will *never* be linear, and its gradient 
never be constant on cells. It would not surprise me if this yielded 
issues in the convergence theory that destroy your convergence order. 
Indeed, I have seen statements like this in the literature, and it can 
be solved by not using a bilinear mapping from the reference cell to the 
real cell but just a linear mapping determined by 3 of the 4 vertices of 
the cell. This would be bad for continuous finite elements because these 
functions would not be continuous along two of the four edges, but that 
really doesn't matter for DG methods :-)

I'm not enough of a theorist to tell you where to look for these kinds 
of statements, but I'm not surprised by the behavior you describe either.

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] Diffusion equation with DG fails with distorted mesh

2019-04-26 Thread Doug
Hello Gary,

Since we have two independent implementation of DG, I am curious to know if 
you achieve optimal error convergence order of (p+1) on *distorted* grids 
for either linear advection, diffusion, or convection-diffusion. 

Best Regards,

Doug

On Friday, April 26, 2019 at 11:48:58 AM UTC-4, Gary Uppal wrote:
>
> Prof. Bangerth,
>
> Thanks for the reply. I get problems on less distorted meshes as well. 
> For example, in the mesh with a hole shown below.  I eventually need 
> multiple holes and an advection term. Actually, the advection seems ok and 
> I mainly get issues when the diffusion is large. If the distorted mesh is 
> the issue, is there a better way to construct the geometry? Is there not a 
> way to get better accuracy even with the distorted mesh? I've tried 
> refining the mesh and lowering the time step, but that doesn't seem to 
> change the results. 
>
> Thank you,
> Gary   
> [image: with_hole.jpg]
>
> On Fri, 26 Apr 2019 at 04:50, Doug > 
> wrote:
>
>> Prof. Bangerth,
>>
>> I am encountering a similar "issue" regarding distorted grids. I am using 
>> DG for linear advection with a manufactured solution to verify the code.
>>
>> I recover optimal convergence orders (p+1) for straight meshes when 
>> refining in h. I recover optimal orders (p+1) for a sine-transformed grid 
>> with some skewed angles up to around 45 degrees (as seen attached).
>>
>> However, I lose my optimal (p+1) order as soon as I apply any small 
>> amount of random distortion to my final refined grid (see attached). Note 
>> that the distortion factor of 0.15 is used, but this loss of order also 
>> occurs for small distortion factors down to 0.01. Instead of (p+1), I may 
>> recover (p) or even (p-1). This only applies when the final grid is 
>> randomly distorted. If I start from a coarse-ish distorted grid (say 100 
>> cells), and refine globally such that its children aren't distorted, I then 
>> recover my optimal orders again.
>>
>> Since my distorted cell still have reasonable angles no matter its 
>> refinement level, degenerate cells wouldn't be an issue. Also, since I am 
>> not performing adaptive refinement, there are no hanging nodes such that 
>> the triangulation becomes irregular.
>>
>> Would you know happen to know if some other condition is violated from 
>> random distortions such that error estimates do not hold?
>>
>> Best regards,
>>
>> Doug
>>
>> On Thursday, April 25, 2019 at 11:15:31 PM UTC-4, Wolfgang Bangerth wrote:
>>>
>>> On 4/25/19 2:50 PM, Gary Uppal wrote: 
>>> > 
>>> > I am trying to solve the diffusion equation with Discontinuous 
>>> Galerkin 
>>> > elements. The solution looks good with a regular structured mesh, but 
>>> if I 
>>> > distort the mesh, the solution blows up and does not converge. Is 
>>> there an 
>>> > obvious reason this would happen? I later need a mesh with holes in 
>>> it, so I 
>>> > cannot always use the structured mesh. 
>>> > 
>>> > I use the interior penalty method and get the diffusion matrix using 
>>> > MeshWorker as in Step-39. I compute the mass matrix and solve the 
>>> diffusion 
>>> > equation with an implicit backward Euler method and am using periodic 
>>> boundary 
>>> > conditions. Snapshots of the solution with structured and distorted 
>>> meshes are 
>>> > shown below. Any help is appreciated! 
>>>
>>> Gary -- the wrong solutions happen on cells that are very nearly 
>>> degenerate 
>>> (i.e., are almost triangular). The finite element theory says that the 
>>> error 
>>> between the exact and the numerical solution is bounded by a constant 
>>> times 
>>> some power of h, where the constant depends on the minimal and maximal 
>>> angles 
>>> at the vertices of the cells. Theory then also predicts that this 
>>> constant 
>>> goes to infinity if the maximal angle at one vertex comes close to 180 
>>> degrees 
>>> -- which is exactly what is happening in your case. 
>>>
>>> So choose a mesh that is less distorted and you should be fine. 
>>>
>>> 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 a topic in the 
>> Google Groups "deal.II User Group" group.
>> To unsubscribe from this topic, visit 
>> https://groups.google.com/d/topic/dealii/ZSHIDbp7yfE/unsubscribe.
>> To unsubscribe from this group and all its topics, send an email to 
>> dea...@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 

Re: [deal.II] Re: Integrating a deal.II-specific function in a NOX-Interface for MPI-threads

2019-04-26 Thread Wolfgang Bangerth
On 4/26/19 1:59 AM, 'Maxi Miller' via deal.II User Group wrote:
> I tried to reduce it as much as possible while still being able to run, 
> the result is attached. It still has ~600 lines, but I am not sure if I 
> should remove things like the output-function? Furthermore the NOX-logic 
> takes quite a lot of lines, and I would like to keep the program as 
> usable as possible while removing unneccessary lines. Is that example 
> usable?

Remove everything you can. If you get an error at one point, no output 
will be created -- so remove the output function. If the error happens 
before the linear system is solved, remove the code that computes the 
entries of the matrix. Then remove the matrix object itself. If you are 
currently solving a linear system before the error happens, try what 
happens if you just comment out the solution step -- and if the error 
still happens (because it probably doesn't depend on the actual values 
of the vector), then remove the solution step and the assembly of the 
matrix and the matrix. If you need a DoFHandler to set up the vector, 
output the index sets and sizes you use for the processors involved and 
build these index sets by hand instead -- then remove the DoFHandler and 
the finite element and everything else related to this. Continue this 
style of work until you really do have a small testcase :-)

Cheers
  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] Imposing the Neumann BC using right_hand_side() function

2019-04-26 Thread Muhammad Mashhood
It seems that "tmp" vector stays unchanged i.e. function is seemingly not 
assigning any face rhs values to "tmp" vector. 

On Friday, April 26, 2019 at 11:42:09 AM UTC+2, Muhammad Mashhood wrote:
>
> Thank you Prof. Bangerth for the correction that was quit helping. 
> Currently the code started running but unfortunately the function is 
> returning with the zero "tmp" vector. Any possible correction? 
>
>
>
>
>
>
> *tmp.reinit (solution.size());
> VectorTools::create_boundary_right_hand_side(dof_handler,  
>  QGauss(fe.degree+1),
> 
>  ConstantFunction(1), 
>  tmp);*
>
> Currently I also did not specify the id of the node where I want to 
> impose  "*ConstantFunction(1)*" value 1 supposing that the "
> *create_boundary_right_hand_side*" will apply it itself on both ends when 
> ids not specified explicitly. 
>
> On Thursday, April 25, 2019 at 8:49:04 PM UTC+2, Wolfgang Bangerth wrote:
>>
>> On 4/25/19 9:46 AM, Muhammad Mashhood wrote: 
>> > 
>> > As unlike the step-8 here in step-26 the rhs is being formed directly 
>> > globally as "system_rhs" using "create_right_hand_side()" function so 
>> > can I use the "create_boundary_right_hand_side()" function using 
>> > following lines of code? 
>>
>> It is definitely possible to *add* to the vector you get from 
>> create_boundary_rhs(). But in this call... 
>>
>> > / 
>> > / 
>> > /std::map::iterator 
>> > it=boundary_count.begin(); 
>> > 
>> >  VectorTools::create_boundary_right_hand_side(dof_handler, 
>> >   
>> >   QGauss(fe.degree+1), 
>> >   rhs_function, 
>> >   tmp,/ 
>> > / it->first);/ 
>> > where the "rhs_function" provides the values of /g0/ and /g1 /at 
>> > boundary node points and "it->first" is the boundary id where I want to 
>> > apply the BC. 
>>
>> ...you are asked to pass in a Quadrature, not a Quadrature. 
>> Consequently, the code does not compile. This is actually stated quite 
>> explicitly in the second part of the error message: 
>>
>> dealii::VectorTools::create_boundary_right_hand_side(const 
>> dealii::DoFHandler&, const dealii::Quadrature<(dim - 
>> 1)>&, const dealii::Function> VectorType::value_type>&, VectorType&, const std::set&) 
>> [with int dim = 1; int spacedim = 1; VectorType = 
>> dealii::Vector; typename VectorType::value_type = double] 
>> void create_boundary_right_hand_side 
>>  ^ 
>> /usr/local/include/deal.II/numerics/vector_tools.h:2105:8: note:   no 
>> known conversion for argument 2 from ‘dealii::QGauss<1>’ to ‘const 
>> dealii::Quadrature<0>&’ 
>>
>>
>> Best 
>>   WB 
>>
>> -- 
>>  
>> 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] Imposing the Neumann BC using right_hand_side() function

2019-04-26 Thread Muhammad Mashhood
Thank you Prof. Bangerth for the correction that was quit helping. 
Currently the code started running but unfortunately the function is 
returning with the zero "tmp" vector. Any possible correction? 






*tmp.reinit (solution.size());
VectorTools::create_boundary_right_hand_side(dof_handler,  
 QGauss(fe.degree+1),

 ConstantFunction(1), 
 tmp);*

Currently I also did not specify the id of the node where I want to impose  
"*ConstantFunction(1)*" value 1 supposing that the "
*create_boundary_right_hand_side*" will apply it itself on both ends when 
ids not specified explicitly. 

On Thursday, April 25, 2019 at 8:49:04 PM UTC+2, Wolfgang Bangerth wrote:
>
> On 4/25/19 9:46 AM, Muhammad Mashhood wrote: 
> > 
> > As unlike the step-8 here in step-26 the rhs is being formed directly 
> > globally as "system_rhs" using "create_right_hand_side()" function so 
> > can I use the "create_boundary_right_hand_side()" function using 
> > following lines of code? 
>
> It is definitely possible to *add* to the vector you get from 
> create_boundary_rhs(). But in this call... 
>
> > / 
> > / 
> > /std::map::iterator 
> > it=boundary_count.begin(); 
> > 
> >  VectorTools::create_boundary_right_hand_side(dof_handler, 
> >   
> >   QGauss(fe.degree+1), 
> >   rhs_function, 
> >   tmp,/ 
> > / it->first);/ 
> > where the "rhs_function" provides the values of /g0/ and /g1 /at 
> > boundary node points and "it->first" is the boundary id where I want to 
> > apply the BC. 
>
> ...you are asked to pass in a Quadrature, not a Quadrature. 
> Consequently, the code does not compile. This is actually stated quite 
> explicitly in the second part of the error message: 
>
> dealii::VectorTools::create_boundary_right_hand_side(const 
> dealii::DoFHandler&, const dealii::Quadrature<(dim - 
> 1)>&, const dealii::Function VectorType::value_type>&, VectorType&, const std::set&) 
> [with int dim = 1; int spacedim = 1; VectorType = 
> dealii::Vector; typename VectorType::value_type = double] 
> void create_boundary_right_hand_side 
>  ^ 
> /usr/local/include/deal.II/numerics/vector_tools.h:2105:8: note:   no 
> known conversion for argument 2 from ‘dealii::QGauss<1>’ to ‘const 
> dealii::Quadrature<0>&’ 
>
>
> Best 
>   WB 
>
> -- 
>  
> 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] Re: Integrating a deal.II-specific function in a NOX-Interface for MPI-threads

2019-04-26 Thread 'Maxi Miller' via deal.II User Group
I tried to reduce it as much as possible while still being able to run, the 
result is attached. It still has ~600 lines, but I am not sure if I should 
remove things like the output-function? Furthermore the NOX-logic takes 
quite a lot of lines, and I would like to keep the program as usable as 
possible while removing unneccessary lines. Is that example usable?

Am Freitag, 26. April 2019 05:12:39 UTC+2 schrieb Wolfgang Bangerth:
>
> On 4/25/19 10:13 AM, 'Maxi Miller' via deal.II User Group wrote: 
> > After running some tests, I ended up reducing the problem to the 
> transfer to 
> > and from the Epetra-vectors. I have to write an interface to the model 
> > (according to the instructions), and as shown in the code above. There I 
> have 
> > Epetra-Vectors created by my interface, with a size of 
> > locally_relevant_dofs.size(), and TrilinosWrappers::MPI::Vectors. 
> Transfer 
> > from the Epetra-Vectors to the MPI::Vectors works fine, but the way back 
> does 
> > not work (The MPI::Vectors are larger than the Epetra_Vectors). Are 
> there any 
> > hints in how I still could fit the data from the MPI::Vector into the 
> > Epetra_Vector? Or should I rather ask this on the Trilinos mailing list? 
>
> Good question. I think it would probably be very useful to have small 
> testcase 
> others could look at. The program you have attached is still 1,300 lines, 
> which is far too much for anyone to understand. But since you have an idea 
> of 
> what the problem is, do you think you could come up with a small testcase 
> that 
> illustrates exactly the issue you have? It doesn't have to do anything 
> useful 
> at all -- in your case, I think all that's necessary is to create a 
> vector, 
> convert it as you describe above, and then output the sizes to show that 
> they 
> are not the same. Run this on two processors, and you should have 
> something 
> that could be done in 50 or 100 lines, and I think that might be small 
> enough 
> for someone who doesn't know your code to understand what is necessary. 
>
> I've built these sorts of testcases either from scratch in the past, or by 
> just keep removing things from a program that (i) are run after the time 
> the 
> problem happens, and (ii) then removing everything that is not necessary. 
>
> 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.
//Nox include files
#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 

// NOX Objects
#include "NOX.H"
#include "NOX_Epetra.H"

// Trilinos Objects
#ifdef HAVE_MPI
#include "Epetra_MpiComm.h"
#else
#include "Epetra_SerialComm.h"
#endif
#include "Epetra_Vector.h"
#include "Epetra_Operator.h"
#include "Epetra_RowMatrix.h"
#include "NOX_Epetra_Interface_Required.H" // base class
#include "NOX_Epetra_Interface_Jacobian.H" // base class
#include "NOX_Epetra_Interface_Preconditioner.H" // base class
#include "Epetra_CrsMatrix.h"
#include "Epetra_Map.h"
#include "Epetra_LinearProblem.h"
#include "AztecOO.h"
#include "Teuchos_StandardCatchMacros.hpp"

// User's application specific files
//#include "problem_interface.h" // Interface file to NOX
//#include "finite_element_problem.h"

// Required for reading and writing parameter lists from xml format
// Configure Trilinos with --enable-teuchos-extended
#ifdef HAVE_TEUCHOS_EXTENDED
#include "Teuchos_XMLParameterListHelpers.hpp"
#include 
#endif

// This is C++ ...
#include 
#include 

using namespace dealii;

template 
class BoundaryValues : public dealii::Function
{
public:
BoundaryValues(const size_t n_components)
: dealii::Function(1), n_components(n_components)
{}

virtual double value (const dealii::Point   ,
  const unsigned int  component = 0) const;

virtual void vector_value(const dealii::Point , dealii::Vector ) const;
private:
const size_t n_components;
};


template 
double BoundaryValues::value (const dealii::Point ,
   const unsigned int) const
{
const double x = p[0];
const double y 

[deal.II] Re: Post-processing with different cell material properties

2019-04-26 Thread Tuanny Cajuhi
Dear all,

can I still use this implementation to output the material_ids together 
with the updated DataPostprocessor class?

Thanks.

Best,
Tuanny 



On Sunday, August 14, 2016 at 3:54:18 PM UTC+2, Pete Griffin wrote:
>
> Thank you Daniel for your response.
>
> I took a while since I was trying to implement something that was similar 
> to the post you pointed to by Philip Wardlaw, but was as requested by 
> Wolfgang's response a little more general. I looked on git but did not know 
> how to search to find what he may or may not have posted by him.
>
> I have minimally modified four files that are attached (data_out.h, 
> data_postprocessor.h, data_out.cc and data_postprocessor.cc). The 
> modifications start and end with "// Start" or "// End". Their is 
> an addition of an enum in data_out.h. The same could be done for the other 
> places build_patches() is used like data_out_faces() etc. but I am not sure 
> whether would be useful since I am not familiar with them.
>
>   typedef enum
>   {
> no_data,
> material_id_data,
> user_index_data,
> pointer_data, // not implemented
>   } ExtraDataOut;
>
>
> The value of the enum is passed from user code to build_patches () as an 
> unsigned int. I guess ideally it would be an DataOut::ExtraDataOut but 
> my C++ was not good enough to have the enum known to 
> DataPostprocessor. The new call to build_patches() would be:
>
> for HP code:
>   data_out.build_patches (ndiv, DataOut::user_index_data); // not yet 
> mapping_collection
> or NON-HP code:
>   data_out.build_patches (mapping, ndiv, DataOut::curved_inner_cells, 
> DataOut::user_index_data);
>
> In both cases the final new parameter defalts to DataOut::no_data. I 
> created two new functions:
>   compute_derived_quantities_scalar_with_data ()
>   compute_derived_quantities_vector_with_data ()
> with the data item returned as an unsigned long long int to accomadate a 
> 64 bit pointer and kept the original ones as is:
>   compute_derived_quantities_scalar ()
>   compute_derived_quantities_vector ()
>
> When attempting to have just one I got the warning related to "hidden 
> funtion" that I was not able to get rid of so I created another name for 
> the new one which for me worked fine.
>
> I have "make clean debug" and "make release run" the following steps:
>   step-18, 17, 8, 27, 30, 6 and 51 and my extended step-8 with HP, NON-HP 
> and MPI/PETSc using the original compute_derived_quantities_vector() and 
> new compute_derived_quantities_vector_with_data() without issues. All 
> extended step-8 gave the same numerical results, considering HP/NON-HP.
>
> I hope this will be useful to others.
>
> Pete Griffin
>
>
> On Saturday, August 6, 2016 at 4:29:46 PM UTC-4, Daniel Arndt wrote:
>>
>> Pete,
>>
>> I remember someone having something like this in mind. This is the link 
>> to that discussion:
>> https://groups.google.com/forum/#!topic/dealii/QKuDndKojug
>>
>> A different way would be to derive from DataOut and overwrite 
>> DataOut::first_cell and DataOut::next_cell by skipping all the cells that 
>> don't have a specific material_id.
>> If you do this for each material_id, you could create an output file for 
>> each material_id and join them afterwards. This would be similar to how you 
>> combine output from different MPI processes using a pvtu file (step-40).
>>
>> 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.