[deal.II] Re: Extra layer around mesh

[Joel Davidsson]

Dear Daniel,

Here is the small example you asked for. If you run it as it is, it will 
work. If you can 1 to 2 in line 107, it crashes. Im run dealii-8.2.1, if 
that is needed.

Best,
Joel

On Tuesday, November 1, 2016 at 11:44:48 AM UTC+1, Daniel Arndt wrote:
>
> Joel,
>
> I have code that runs fine with periodic boundary condition for 
>> a GridGenerator::parallelepiped. Also when I change the mesh to 
>> a GridGenerator::subdivided_parallelepiped and has only 1 as 
>> n_subdivisions, it also works. But when I change it to 2, I get the 
>> following error:
>> [...]
>>
> GridTools::collect_periodic_faces[1] tries to match the faces on the 
> boundaries specified by the two boundary_ids on the coarsest level and here 
> finds a face that can't be matched. 
> Can you give us a minimal example just creating that mesh and trying to 
> compute periodic boundary conditions?
> From just your description, it is hard to say what the problem is, but I 
> would not expect that subdividing should be an issue.
>
> Best,
> Daniel
>
> [1] 
> https://dealii.org/8.4.1/doxygen/deal.II/namespaceGridTools.html#aaeadfc0053429f542fbfd48d192b94f0
>

-- 
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) 1999 - 2013 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 at
 * the top level of the deal.II distribution.
 *
 * -

 *
 * Author: Wolfgang Bangerth, University of Heidelberg, 1999
 */

#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 Step4
{
public:
  Step4 ();
  void run ();

private:
  void make_grid ();
  void setup_system();
  void assemble_system ();
  void solve ();
  //void output_results () const;

  Triangulation   triangulation;
  FE_Qfe;

  DoFHandler  dof_handler;

  SparsityPattern  sparsity_pattern;
  SparseMatrix system_matrix;

  ConstraintMatrix constraints;

  Vector   solution;
  Vector   system_rhs;
};

// basis
Point<3> basis1(4.0,0.0,0.0);
Point<3> basis2(0.0,4.0,0.0);
Point<3> basis3(0.0,0.0,4.0);

template 
Step4::Step4 ()
  :
  fe (1),
  dof_handler (triangulation)
{}

template 
void Step4::make_grid ()
{
	Point basis [dim] = {basis1,basis2,basis3};

	// change 1 to 2 -> fail
	GridGenerator::subdivided_parallelepiped(triangulation,1,basis,true);
	//GridGenerator::parallelepiped(triangulation,basis,true);

	// Translate the main grid so that midpoint is the center
	const Point translation = 0.5*(basis1+basis2+basis3);
	GridTools::shift(-translation,triangulation);

	// refinement
	//triangulation.refine_global(2);

	std::ofstream out ("./grid.vtk");
	GridOut grid_out;
	grid_out.write_vtk (triangulation, out);
	std::cout << "Grid has been written to grid.vtk" << std::endl;
}

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

	std::cout << "   Number of degrees of freedom: "
	  << dof_handler.n_dofs()
	  << std::endl;

	constraints.clear ();

	bool periodic = true;
	bool x_dir = true;
	bool y_dir = true;
	bool z_dir = true;

	if (periodic == true)
	{
		if (x_dir == true)
		{
			Tensor<1,dim> offset;
			offset[0]=basis1[0];
			offset[1]=basis1[1];
			offset[2]=basis1[2];
			DoFTools::make_periodicity_constraints(dof_handler,0,1,0,offset,constraints);
		}
		else if (x_dir == false)
		{
			VectorTools::interpolate_boundary_values (dof_handler,0,ZeroFunction(),constraints);
			VectorTools::interpolate_boundary_values (dof_handler,1,ZeroFunction(),constraints);
		}

		if (y_dir == true)
		{
			Tensor<1,dim> offset;
			offset[0]=basis2[0];
			offset[1]=basis2[1];
			offset[2]=basis2[2];
			DoFTools::make_periodicity_constraints(dof_handler,2,3,1,offset,constraints);
		}
		else if (y_dir == false)
		{
			VectorTools::interpolate_boundary_values (dof_handler,2,ZeroFunction(),constraints);
			VectorTools::interpolate_boundary_values 

[deal.II] Re: Extra layer around mesh

[Joel Davidsson]

Dear all,

I have code that runs fine with periodic boundary condition for 
a GridGenerator::parallelepiped. Also when I change the mesh to 
a GridGenerator::subdivided_parallelepiped and has only 1 as 
n_subdivisions, it also works. But when I change it to 2, I get the 
following error:

*An error occurred in line <2686> of file 
 in function*
*void dealii::GridTools::collect_periodic_faces(const CONTAINER&, 
dealii::types::boundary_id, dealii::types::boundary_id, int, 
std::vector<dealii::GridTools::PeriodicFacePair >&, const dealii::Tensor<1, CONTAINER:: 
space_dimension>&, const dealii::FullMatrix&, const 
std::vector&) [with CONTAINER = dealii::DoFHandler<3>; 
dealii::types::boundary_id = unsigned char; typename 
CONTAINER::cell_iterator = 
dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<3>, false> 
>]*
*The violated condition was: *
*pairs1.size() == pairs2.size()*
*The name and call sequence of the exception was:*
*ExcMessage ("Unmatched faces on periodic boundaries")*
*Additional Information: *
*Unmatched faces on periodic boundaries*

How can I fix it?

Best,
Joel

On Friday, August 12, 2016 at 9:35:23 AM UTC+2, Joel Davidsson wrote:
>
> Ok
>
> Thanks,
>
> Joel
>
> On Thursday, August 11, 2016 at 6:05:41 PM UTC+2, Bruno Turcksin wrote:
>>
>> Joel,
>>
>> On Thursday, August 11, 2016 at 11:29:35 AM UTC-4, Joel Davidsson wrote:
>>>
>>>
>>> I just had a question about the third option. If one used FE_Nothing, 
>>> what is the boundary that is used when doing calculations? Is it the 
>>> boundary of the big box or the green box? Is this handle automatically or 
>>> does one need to mark the green box with boundary indicator for example.
>>>
>> It will be the boundary on the big box and you can't use boundary 
>> indicator because you are not on the boundary. The boundary is a geometric 
>> "quantity", it doesn't care about the type of finite elements you are 
>> using. So if you want to use Dirichlet "boundary" condition on the green 
>> box, you need to build the ConstraintMatrix yourself.
>>
>> Best,
>>
>> Bruno
>>
>

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

[deal.II] Re: Vector-valued gradient of solution vector

[Joel Davidsson]

Dear Daniel,

I did try it, but I cant just change the local_velocity_values from a 
std::vector> to std::vector. Then the 
fe_values[velocities].get_function_values function doesnt compile. So I 
commented that part out, the results you see is from the code as it is, I 
used the for loop at line 354 to get the size (which I know is wrong, but 
at least it give something).  Even if I had found the correct way to do 
this, it would not have solved my original problem.

(correct me if wrong) This method of using norm_sqr() would only produce 
the size, not the individual x,y and z components. I need the x,y and 
z components separately so I can get the size but also the difference 
between two solutions and that error of these two solutions. Hence my 
previous message. I want to be able to manipulate the vector components. 
How would I do this?

Best,
Joel


On Tuesday, September 27, 2016 at 7:38:48 PM UTC+2, Daniel Arndt wrote:
>
> Joel,
>
> Its been a week since my latest post, I was just wondering if you could 
>> help me.
>>
> Did you try what I suggested? What were the problems?
>
> Looking at the latest files you sent, local_velocity_values has still not 
> the type std::vector
> and you are still using
>
> for (unsigned int q_index=0; q_index   for (unsigned int i=0; i velocity[local_dof_indices[i]] = local_velocity_values[q_index];
>
> instead of
>
> for (unsigned int i=0; i   velocity[local_dof[indices[i]] = local_velocity_values(i).norm_sqr();
>
> This way you don't rely on any particular order of the dofs.
>
> 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: Vector-valued gradient of solution vector

[Joel Davidsson]

Dear Daniel,

To clarify what I am searching for, I would like to get the vector-valued 
solution into a format that matches a scalar. Either as a 
std::vector> or a std::vector, the format doesn't 
matter as long as it matches a scalar dofs. So I can manipulate the vector 
components and for example get the size.

I made a test program that I have attached, that takes the gradient of a 
gaussian function. At line 354, I do the convergence from a vector valued 
solution to get the size of the vector (as I posted before). This was wrong 
and you guys said that I shouldn't assume any particular order of the dofs. 
Then I tried to this conversion the correct way, see line 365, but I do not 
know how to do this conversion the correct way.

The code I have sent you, produces the correct result for fe1 and fe2, i.e. 
when I set fe_scalar (2), fe(FE_Q(2),dim) at line 100. But failed for 
fe3. I have attached the result I got for fe2 and fe3 as well as the plot 
script I have used to generate these plots.

I would be really grateful if you could help me figure out how to solve 
this.

Best,
Joel

On Monday, September 19, 2016 at 11:56:51 AM UTC+2, Daniel Arndt wrote:
>
> Joel,
>
>> [...]
>>  for (unsigned int q_index=0; q_index>for (unsigned int i=0; i>  {
>> velocity[local_dof_indices[i]] = local_velocity_values[q_index];
>>  }
>>
>> Here you assign the same Tensor<1,dim> at all DoFs multiple times and you 
> end with velocity[*]=local_velocity_values[n_q_points]-1.
> I don't think that this is what you want to do, is it? If I understand you 
> correctly, you try to assign the magnitude of the velocity in each degree 
> of freedom
> to a different scalar Vector that represents a scalar FiniteElement Vector.
> In this case, you need to use a Quadrature object that uses the 
> unit_support_points of the scalar FE [1] and your assignment would read
>
> for (unsigned int i=0; i   velocity[local_dof[indices[i]] = local_velocity_values(i).norm_sqr();
>
> where velocity is a Vector that is related to fe_scalar and 
> local_velocity_values has the type  std::vector.
>
> Best,
> Daniel
>
> [1] 
> https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#how-to-get-the-mapped-position-of-support-points-of-my-element
>

-- 
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) 1999 - 2013 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 at
 * the top level of the deal.II distribution.
 *
 * -

 *
 * Author: Wolfgang Bangerth, University of Heidelberg, 1999
 */

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

const double pi=3.14159265359;
unsigned int center_id=0;

template 
class Step4
{
public:
  Step4 ();
  void run ();

private:
  void make_grid ();
  void setup_system();
  void assemble_system ();
  void solve ();
  void output_results () const;

  Triangulation   triangulation;
  FE_Qfe_scalar;
  FESystem		   fe;

  DoFHandler  dof_handler_scalar;
  DoFHandler  dof_handler;

  SparsityPattern  sparsity_pattern;
  SparseMatrix system_matrix;

  ConstraintMatrix constraints;

  Vector   solution;
  Vector   system_rhs;
  Vector   test;
  Vector   size;

  std::vector nodes;
};

template 
Step4::Step4 ()
  :
  fe_scalar (3),
  fe (FE_Q(3),dim),
  dof_handler (triangulation),
  dof_handler_scalar (triangulation)
{}

template 
void Step4::make_grid ()
{
  GridGenerator::hyper_cube (triangulation, -2.5, 2.5);
  triangulation.refine_global (2);

  for (unsigned int step=0; step<2; ++step)
  {
	  typename Triangulation::active_cell_iterator
	  cell = triangulation.begin_active(),
	  endc = triangulation.end();
	  for (; cell!=endc; ++cell)
	  {
		  for (unsigned int v=0;v < 

[deal.II] Re: Vector-valued gradient of solution vector

[Joel Davidsson]

Dear J-P,

Ok, I tried to to use the FEValueExtractors but I am not sure if my code is 
working correctly. 

   std::vector<Tensor<1,dim>>   velocity;
   velocity.resize (dof_handler_scalar.n_dofs());

 QGauss  quadrature_formula(dof_handler.get_fe().degree+1);

 FEValues fe_values (fe, quadrature_formula,
  update_values   | update_gradients |
  update_quadrature_points | update_JxW_values);

 FEValues fe_values_scalar (fe_scalar, quadrature_formula,
  update_values   | update_gradients |
  update_quadrature_points | update_JxW_values);

 const unsigned int   dofs_per_cell = fe_scalar.dofs_per_cell;
 const unsigned int   n_q_points= quadrature_formula.size();

 std::vector local_dof_indices (dofs_per_cell);
  std::vector<Tensor<1,dim> > local_velocity_values (n_q_points);

 const FEValuesExtractors::Vector velocities (0);

 typename DoFHandler::active_cell_iterator
   cell_scalar = dof_handler_scalar.begin_active(),
   endc_scalar = dof_handler_scalar.end();

 typename DoFHandler::active_cell_iterator
 cell = dof_handler.begin_active(),
 endc = dof_handler.end();

 for (; cell!=endc; ++cell, ++cell_scalar)
   {
 fe_values.reinit (cell);
 fe_values_scalar.reinit (cell_scalar);
 fe_values[velocities].get_function_values 
(solution,local_velocity_values);

 cell_scalar->get_dof_indices (local_dof_indices);

 for (unsigned int q_index=0; q_index<n_q_points; ++q_index)
   for (unsigned int i=0; i<dofs_per_cell; ++i)
 {
velocity[local_dof_indices[i]] = local_velocity_values[q_index];
 }
 }

Best,
Joel

On Tuesday, September 13, 2016 at 2:19:43 PM UTC+2, Jean-Paul Pelteret 
wrote:
>
> Dear Joel,
>
> In general, no you cannot. Here you assume a specific ordering of the DoFs 
> in the global solution vector, and this changes for each finite element. 
> For example, using FE_Q's might result the ordering [x0 y0 z0 x1 y1 z1 ...] 
> and FE_DGQ's might lead to [x0 x1 x2 ... y0 y1 y2 ... z0 z1 z2...]. To do 
> this safely you should definitely do this using the built in tools, such as 
> FEValueExtractors 
> <https://www.dealii.org/developer/doxygen/deal.II/namespaceFEValuesExtractors.html>
>  
> or  fe.system_to_component_index 
> <https://www.dealii.org/developer/doxygen/deal.II/classFiniteElement.html#a27220a135402b96c7e6eecbb04acda56>.
>  
> See the tutorials for examples on how to use them.
>
> Regards,
> J-P
>
> On Tuesday, September 13, 2016 at 2:05:56 PM UTC+2, Joel Davidsson wrote:
>>
>> Dear all,
>>
>> If I have a Vector-valued solution in 3D and I would like to get the size 
>> of each vector in another smaller Vector, can I do it like this:
>>
>> for (unsigned int i=0;i<solution.size();i=i+3)
>> {
>>  unsigned int x=i;
>>  unsigned int y=i+1;
>>  unsigned int z=i+2;
>>  unsigned int j = i / 3.0;
>>
>>  size(j)=pow(solution(x),2.0)+pow(solution(y),2.0)+pow(solution(z),2.0);
>> }
>>
>> Where solution is the vector valued Vector that is 3 times bigger that 
>> the scalar Vector called size. 
>>
>> Best,
>> Joel
>>
>> On Tuesday, August 30, 2016 at 6:48:46 PM UTC+2, Jean-Paul Pelteret wrote:
>>>
>>> Dear Joel,
>>>
>>> No, this functionality doesn't exist. However, you can create a 
>>> quadrature rule with the quadrature points located at the support points 
>>> (see the FAQ 
>>> <https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#how-to-get-the-mapped-position-of-support-points-of-my-element>),
>>>  
>>> and could then extract the solution based on your custom quadrature rule.
>>>
>>> Regards,
>>> J-P
>>>
>>> On Tuesday, August 30, 2016 at 5:32:13 PM UTC+2, Joel Davidsson wrote:
>>>>
>>>> Dear Daniel,
>>>>
>>>> I was just wondering if their exist a similar function 
>>>> to DoFTools::map_dofs_to_support_points, but that maps a vector-valued 
>>>> solution Vector into to a std::vector<Point> instead?
>>>>
>>>> Best,
>>>> Joel
>>>>
>>>> On Monday, August 22, 2016 at 9:03:22 PM UTC+2, Daniel Arndt wrote:
>>>>>
>>>>> Joel,
>>>>>
>>>>> You should not rely on any particular sorting of the dofs in the 
>>>>> solution vector. Instead you can ask the FiniteElement on each cell for 
>>>>> the 
>>>>> component of a local dof by
>>>>> const unsigned int component = fe.system_to_component_index(i).first; 
>>

[deal.II] Re: Vector-valued gradient of solution vector

[Joel Davidsson]

Dear all,

If I have a Vector-valued solution in 3D and I would like to get the size 
of each vector in another smaller Vector, can I do it like this:

for (unsigned int i=0;i<solution.size();i=i+3)
{
 unsigned int x=i;
 unsigned int y=i+1;
 unsigned int z=i+2;
 unsigned int j = i / 3.0;

 size(j)=pow(solution(x),2.0)+pow(solution(y),2.0)+pow(solution(z),2.0);
}

Where solution is the vector valued Vector that is 3 times bigger that the 
scalar Vector called size. 

Best,
Joel

On Tuesday, August 30, 2016 at 6:48:46 PM UTC+2, Jean-Paul Pelteret wrote:
>
> Dear Joel,
>
> No, this functionality doesn't exist. However, you can create a quadrature 
> rule with the quadrature points located at the support points (see the FAQ 
> <https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#how-to-get-the-mapped-position-of-support-points-of-my-element>),
>  
> and could then extract the solution based on your custom quadrature rule.
>
> Regards,
> J-P
>
> On Tuesday, August 30, 2016 at 5:32:13 PM UTC+2, Joel Davidsson wrote:
>>
>> Dear Daniel,
>>
>> I was just wondering if their exist a similar function 
>> to DoFTools::map_dofs_to_support_points, but that maps a vector-valued 
>> solution Vector into to a std::vector<Point> instead?
>>
>> Best,
>> Joel
>>
>> On Monday, August 22, 2016 at 9:03:22 PM UTC+2, Daniel Arndt wrote:
>>>
>>> Joel,
>>>
>>> You should not rely on any particular sorting of the dofs in the 
>>> solution vector. Instead you can ask the FiniteElement on each cell for the 
>>> component of a local dof by
>>> const unsigned int component = fe.system_to_component_index(i).first; 
>>>
>>> Apart from that the DataOut object knows about components and you can do 
>>> some postprocessing on the solution vector using DataPostprocessor [1]
>>>
>>> Best,
>>> Daniel
>>>
>>>  [1] 
>>> https://dealii.org/8.4.1/doxygen/deal.II/classDataPostprocessor.html
>>>
>>> Am Montag, 22. August 2016 17:05:08 UTC+2 schrieb Joel Davidsson:
>>>>
>>>> Dear Daniel,
>>>>
>>>> Thank you for a very good answer, adding the fe_values_scalar and 
>>>> cell_scalar fixed the problem.
>>>>
>>>> I have a follow-up question about the solution I get out, how is the 
>>>> data organized in the solution vector? Say for example I want to loop over 
>>>> all the x components, how would I do that?
>>>>
>>>> Best,
>>>> Joel
>>>>
>>>> On Friday, August 12, 2016 at 5:35:07 PM UTC+2, Daniel Arndt wrote:
>>>>>
>>>>> Joel,
>>>>>
>>>>> Yes the matrix you are assembling is a vector-valued mass matrix now.
>>>>> For me your code is failing with
>>>>>
>>>>> void dealii::FEValuesBase<dim, spacedim>::get_function_gradients(const 
>>>>> InputVector&, std::vector<dealii::Tensor<1, spacedim, typename 
>>>>> InputVector::value_type> >&) const [with InputVector = 
>>>>> dealii::Vector; int dim = 3; int spacedim = 3; typename 
>>>>> InputVector::value_type = double]
>>>>> The violated condition was: 
>>>>> (fe->n_components()) == (1)
>>>>> The name and call sequence of the exception was:
>>>>> dealii::ExcDimensionMismatch((fe->n_components()),(1))
>>>>> Additional Information: 
>>>>> Dimension 3 not equal to 1
>>>>>
>>>>> and this is not surprising. What you are doing is to extract the 
>>>>> gradients of a vector-valued finite element solution. The object that 
>>>>> should store this should therefore be a 
>>>>> std::vector<std::vector<Tensor<1,dim>> as you want to store for each 
>>>>> quadrature point and each component a Tensor<1,dim>.
>>>>>
>>>>> What you want to do is really that in "Do not work". As you have two 
>>>>> DoFHandlers, you should also have two FEValues objects and two 
>>>>> cell_iterator corresponding to the correct DoFHandler. Then you can 
>>>>> extract 
>>>>> the gradient of your scalar field and project it onto the ansatz space 
>>>>> for 
>>>>> the vector-valued field. This should look like:
>>>>>
>>>>> typename DoFHandler::active_cell_iterator cell_scalar = 
>>>>> dof_handler_scalar.begin_active();
>>>>> typename DoFHandler::active_cell_iterator cell = 
>>>>> dof_handler.begin_active();
>>>>> ...
>>>>>
>>>>>  for (; cell!=endc; ++cell, ++cell_vector)
>>>>> {
>>>>>   fe_values.reinit (cell);
>>>>>   fe_values_scalar.reinit (cell_scalar);
>>>>>   fe_values_scalar.get_function_gradients(test,fe_function_grad);
>>>>>   ...
>>>>> }
>>>>>
>>>>> 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: nearest neighbor

[Joel Davidsson]

Thanks Daniel,

I had confused the dof_index with the vertex_number. Now It works as 
expected.

Best,
Joel

On Wednesday, September 7, 2016 at 3:47:21 PM UTC+2, Daniel Arndt wrote:
>
> Joel,
>
> In GridTools::find_cells_adjacent_to_vertex(dof_handler,center_id) you 
> need to give the number of the vertex_number as center_id, not the number 
> of the degree of freedom.
> Then, you want to use 
>
> Quadrature  
> quadrature_formula(dof_handler.get_fe().get_unit_support_points());
>
> to create a Quadrature object that gives you the support points on each 
> cell.
> Finally, you can output the global dof number and their support point on 
> each cell by: 
>
> for (unsigned int j=0; j   std::cout << local_dof_indices[j] << ": "<< 
> fe_values.quadrature_point(j) << std::endl;
>
> In this ouput, you can then search for the dofs with support_point nearest 
> to the vertex you wanted to consider.
>
> 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] nearest neighbor

[Joel Davidsson]

Dear all,

If I have the index to a dof, is there any easy way I can get indexes to 
the nearest neighbor?

Thanks,

Joel

-- 
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: Vector-valued gradient of solution vector

[Joel Davidsson]

Dear Daniel,

Thank you for a fast answer. I think I got the matrix and right-hand side 
correct now, please see test.cc for the example code.

What I want to from the start was to get the vector-valued gradient from a 
scalar field, but I could not get this to work. See section "Do not work" 
in function setup_system in test.cc, this gives 0 everywhere. I want to use 
the information from test vector but since it was made with a different 
dof_handler it doesn´t seem to work. Is there a way around this problem? 
Can a vector be transferred from one dof_handler to another? or is there a 
better way to solve this?

Best,
Joel

On Friday, August 12, 2016 at 12:00:36 AM UTC+2, Daniel Arndt wrote:
>
> Joel,
>
> Does the equation you want to solve for have multiple components? 
> Otherwise it would not make sense to multiply in the right hand side 
> something with a gradient.
> If you want to project the gradient into a vector valued finite element 
> ansatz space, then the matrix you are assembling looks weird.
> In what you are doing you are also considering the coupling between all 
> components.
>
> You are calculating the gradient of your scalar finite element function 
> correctly, but you need to find the correct component to use.
> For doing this you can use something like:
>
> const unsigned int component = fe.system_to_component_index(i).first;
> cell_rhs(i) += ((fe_values.shape_value (i, q_index) *
> fe_function_grad[q_index][component])*
> fe_values.JxW (q_index));
>
> If you want to assemble a mass matrix for a vector-valued finite element, 
> you have also to restrict assembling for the matrix to the case
> where the component for the ith and jth ansatz function are the same.
>
> You might want to have a look at step-8 [1] and the module "Handling 
> vector valued problems"[2].
>
> Best,
> Daniel
>
> [1] 
> https://www.dealii.org/8.4.1/doxygen/deal.II/step_8.html#ElasticProblemassemble_system
> [2] 
> https://www.dealii.org/8.4.1/doxygen/deal.II/group__vector__valued.html
>

-- 
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) 1999 - 2013 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 at
 * the top level of the deal.II distribution.
 *
 * -

 *
 * Author: Wolfgang Bangerth, University of Heidelberg, 1999
 */

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

const double pi=3.14159265359;
double center_id=0;

template 
class Step4
{
public:
  Step4 ();
  void run ();

private:
  void make_grid ();
  void setup_system();
  void assemble_system ();
  void solve ();
  void output_results () const;

  Triangulation   triangulation;
  FE_Qfe_scalar;
  FESystem		   fe;

  DoFHandler  dof_handler_scalar;
  DoFHandler  dof_handler;

  SparsityPattern  sparsity_pattern;
  SparseMatrix system_matrix;

  ConstraintMatrix constraints;

  Vector   solution;
  Vector   system_rhs;
  Vector   test;

  std::vector nodes;
};

template 
Step4::Step4 ()
  :
  fe_scalar (1),
  fe (FE_Q(1),dim),
  dof_handler (triangulation),
  dof_handler_scalar (triangulation)
{}

template 
void Step4::make_grid ()
{
  GridGenerator::hyper_cube (triangulation, -2.5, 2.5);
  triangulation.refine_global (4);

  for (unsigned int step=0; step<2; ++step)
  {
	  typename Triangulation::active_cell_iterator
	  cell = triangulation.begin_active(),
	  endc = triangulation.end();
	  for (; cell!=endc; ++cell)
	  {
		  for (unsigned int v=0;v < GeometryInfo::vertices_per_cell;++v)
		  {
			  Point center;
			  for (unsigned int i=0; ivertex(v));
			  if (std::fabs(distance_from_center) < 1.0)
			  {
  cell->set_refine_flag ();
  break;
			  }
		  }
	  }
	  triangulation.execute_coarsening_and_refinement ();
  }

  std::cout << "   Number of active 

[deal.II] Re: Extra layer around mesh

[Joel Davidsson]

Thank you for a great answer,

The subdivided grid generator approach worked fine for me. 

I just had a question about the third option. If one used FE_Nothing, what 
is the boundary that is used when doing calculations? Is it the boundary of 
the big box or the green box? Is this handle automatically or does one need 
to mark the green box with boundary indicator for example.

Thanks,
Joel

On Thursday, August 4, 2016 at 9:14:18 PM UTC+2, Jean-Paul Pelteret wrote:
>
> Dear Joel,
>
> The way I interpret what you've written here, you're asking two different 
> questions. In the first statement, it sounds like you have a mesh that you 
> wish to increase the size of at some arbitrary time, and the second you 
> know the final bounds *a priori*.
>
> Let me start with your second viewpoint, which is easiest. You could use 
> GridGenerator::subdivided_hyper_cube 
> <https://www.dealii.org/developer/doxygen/deal.II/namespaceGridGenerator.html#adc5d7022d456db0356f11427473f4f76>
>  or GridGenerator::subdivided_hyper_rectangle 
> <https://www.dealii.org/developer/doxygen/deal.II/namespaceGridGenerator.html#ac76417d7404b75cf53c732f456e6e971>
>  to 
> specify any regular discretisation of a Cartesian-aligned volume. So, yes, 
> given that you know the edge length and the size of the elements, you could 
> compute the necessary number of elements in each dimension. If this doesn't 
> give you the exact dimensions of the green box, then you can just move 
> the vertices around 
> <https://www.dealii.org/developer/doxygen/deal.II/step_49.html#grid_3Movingvertices>
>  
> as necessary.
>
> Achieving your goal from the first viewpoint requires more work. To 
> "increase" the size of the grid entails joining triangulations of *zero 
> refinement 
> levels* (i.e. no mesh refinement). So you would need to first create a 
> coarse grid (1 for the green box (which I interpret as being something that 
> you have initially) and then 4 more that share the edges and another 4 that 
> share its vertices, then join them 
> <https://www.dealii.org/developer/doxygen/deal.II/step_49.html#grid_2Mergingtriangulations>.
>  
> Since you have to join equally coarse meshes, all of these meshes probably 
> need to created as "discretised" GridGenerator::subdivided_hyper_rectangle 
> <https://www.dealii.org/developer/doxygen/deal.II/namespaceGridGenerator.html#ac76417d7404b75cf53c732f456e6e971>s
>  
> in the first place.
>
> A third option would be to use a hp::DoFHandler 
> <https://www.dealii.org/developer/doxygen/deal.II/step_46.html> and use 
> FE_Nothing 
> <https://www.dealii.org/developer/doxygen/deal.II/classFE__Nothing.html> 
> elements in the region outside the green box. You could then "grow" the 
> mesh by activating these elements at any time.
>
> I hope that this has been useful to you. Apologies in advance if I don't 
> respond to any further messages - from tomorrow I'll be away from my PC for 
> a while.
>
> Regards,
> J-P
>
> On Thursday, August 4, 2016 at 3:33:29 PM UTC+2, Joel Davidsson wrote:
>>
>> Dear all,
>>
>> Is there some way that one can refine the mesh in the following way, see 
>> attached file: I looking for a way to increase the size of the grid but 
>> keeping the original nodes fixed. If you look in the attached file, you see 
>> that an extra "layer" has formed around the original mesh marked in green. 
>> In this way, the size of the individual element does not change.
>>
>> Another way of looking at it. If I specify a volume and a size of the 
>> smallest element, can I fill this volume with these elements?
>>
>> Best,
>> Joel
>>
>

-- 
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] Extra layer around mesh

[Joel Davidsson]

Dear all,

Is there some way that one can refine the mesh in the following way, see 
attached file: I looking for a way to increase the size of the grid but 
keeping the original nodes fixed. If you look in the attached file, you see 
that an extra "layer" has formed around the original mesh marked in green. 
In this way, the size of the individual element does not change.

Another way of looking at it. If I specify a volume and a size of the 
smallest element, can I fill this volume with these elements?

Best,
Joel

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


refine.pdf
Description: Adobe PDF document