[deal.II] Problem with the installation of deal.II v9.2.0 using Candi

[Hamed Babaei]

Hello everyone,

I'm trying to install the last version of deal.II using candi shell script. 
I have used it multiple times before and it worked well as I install all 
prerequisite packages as it suggests. However, I get the following error 
this time. I appreciate it if anyone could let me know what may be the 
problem.

Thanks and regards,

/home/hbabaei/deal.ii-candi/tmp/unpack/deal.II-v9.2.0/source/dofs/dof_tools_sparsity.cc:1352:55:
 
*error: no matching function for call to ‘make_flux_sparsity_pattern(const 
dealii::hp::DoFHandler<3, 3>&, 
dealii::TrilinosWrappers::BlockSparsityPattern&, 
dealii::AffineConstraints&, const bool&, const dealii::Table<2, 
dealii::DoFTools::Coupling>&, const dealii::Table<2, 
dealii::DoFTools::Coupling>&, const subdomain_id&, )’*

/home/hbabaei/deal.ii-candi/tmp/unpack/deal.II-v9.2.0/source/dofs/dof_tools_sparsity.cc:1352:55:
 
note: candidates are:

/home/hbabaei/deal.ii-candi/tmp/unpack/deal.II-v9.2.0/source/dofs/dof_tools_sparsity.cc:689:3:
 
note: template void 
dealii::DoFTools::make_flux_sparsity_pattern(const DoFHandlerType&, 
SparsityPatternType&)

   make_flux_sparsity_pattern(const DoFHandlerType ,

   ^

/home/hbabaei/deal.ii-candi/tmp/unpack/deal.II-v9.2.0/source/dofs/dof_tools_sparsity.cc:689:3:
 
note:   template argument deduction/substitution failed:

/home/hbabaei/deal.ii-candi/tmp/unpack/deal.II-v9.2.0/source/dofs/dof_tools_sparsity.cc:1352:55:
 
note:   candidate expects 2 arguments, 8 provided

   internal::always_couple_on_faces);

   ^

/home/hbabaei/deal.ii-candi/tmp/unpack/deal.II-v9.2.0/source/dofs/dof_tools_sparsity.cc:520:3:
 
note: template void dealii::DoFTools::make_flux_sparsity_pattern(const 
DoFHandlerType&, SparsityPatternType&, const 
dealii::AffineConstraints&, bool, dealii::types::subdomain_id)

   make_flux_sparsity_pattern(const DoFHandlerType &   dof,

   ^

/home/hbabaei/deal.ii-candi/tmp/unpack/deal.II-v9.2.0/source/dofs/dof_tools_sparsity.cc:520:3:
 
note:   template argument deduction/substitution failed:

/home/hbabaei/deal.ii-candi/tmp/unpack/deal.II-v9.2.0/source/dofs/dof_tools_sparsity.cc:1352:55:
 
note:   cannot convert ‘int_mask’ (type ‘const dealii::Table<2, 
dealii::DoFTools::Coupling>’) to type ‘dealii::types::subdomain_id {aka 
unsigned int}’

   internal::always_couple_on_faces);

   ^

/home/hbabaei/deal.ii-candi/tmp/unpack/deal.II-v9.2.0/source/dofs/dof_tools_sparsity.cc:1334:3:
 
note: template void 
dealii::DoFTools::make_flux_sparsity_pattern(const DoFHandlerType&, 
SparsityPatternType&, const dealii::Table<2, dealii::DoFTools::Coupling>&, 
const dealii::Table<2, dealii::DoFTools::Coupling>&, 
dealii::types::subdomain_id)

   make_flux_sparsity_pattern(const DoFHandlerType &dof,

   ^

/home/hbabaei/deal.ii-candi/tmp/unpack/deal.II-v9.2.0/source/dofs/dof_tools_sparsity.cc:1334:3:
 
note:   template argument deduction/substitution failed:

/home/hbabaei/deal.ii-candi/tmp/unpack/deal.II-v9.2.0/source/dofs/dof_tools_sparsity.cc:1352:55:
 
note:   cannot convert ‘dummy’ (type ‘dealii::AffineConstraints’) 
to type ‘const dealii::Table<2, dealii::DoFTools::Coupling>&’

   internal::always_couple_on_faces);

   ^

/home/hbabaei/deal.ii-candi/tmp/unpack/deal.II-v9.2.0/source/dofs/dof_tools_sparsity.cc:1361:3:
 
note: template void dealii::DoFTools::make_flux_sparsity_pattern(const 
DoFHandlerType&, SparsityPatternType&, const 
dealii::AffineConstraints&, bool, const dealii::Table<2, 
dealii::DoFTools::Coupling>&, const dealii::Table<2, 
dealii::DoFTools::Coupling>&, dealii::types::subdomain_id, const 
std::function&)

   make_flux_sparsity_pattern(

   ^

/home/hbabaei/deal.ii-candi/tmp/unpack/deal.II-v9.2.0/source/dofs/dof_tools_sparsity.cc:1361:3:
 
note:   template argument deduction/substitution failed:

/home/hbabaei/deal.ii-candi/tmp/unpack/deal.II-v9.2.0/source/dofs/dof_tools_sparsity.cc:1352:55:
 
note:   mismatched types ‘const std::function’ and ‘bool(const 
active_cell_iterator&, unsigned int) {aka bool(const 
dealii::TriaActiveIterator, false> >&, unsigned int)}’

   internal::always_couple_on_faces);

   ^

/home/hbabaei/deal.ii-candi/tmp/unpack/deal.II-v9.2.0/source/dofs/dof_tools_sparsity.cc:1352:55:
 
note:   could not resolve address from overloaded function 
‘always_couple_on_faces >’

make[2]: *** 
[source/dofs/CMakeFiles/obj_dofs_release.dir/dof_tools_sparsity.cc.o] Error 
1

make[1]: *** [source/dofs/CMakeFiles/obj_dofs_release.dir/all] Error 2

make: *** [all] Error 2

Failure with exit status: 2

Exit message: There was a problem building dealii v9.2.0.

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

Re: [deal.II] SymmetricTensor rank 6 implementation issue

[Hamed Babaei]

Dear Wolfgang,

I believe it will be less time-consuming if I simply avoid using the 
SymmetricTensor class and use the Tensor class instead and insert the 
symmetric elements manually. 

Thanks anyway for your incredible support!


-- 
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: SymmetricTensor rank 6 implementation issue

[Hamed Babaei]

Dear Bruno,

Thank you very much for your hint. I added the extra include but it did not 
help and I am still getting the same error!

I am wondering what other things I may want to try?

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.

[deal.II] Re: SymmetricTensor rank 6 implementation issue

[Hamed Babaei]

Dear Bruno,

Thank you very much for your hint. I added the extra include but it did not 
help and I am still getting the same error!

I am wondering what other things I may want to try?

Thanks

On Friday, May 3, 2019 at 3:24:04 PM UTC-5, bruno@gmail.com wrote:
>
> Hi,
>
> Try to add #include  We cannot 
> do explicit instantiation for every thing so if you have something that is 
> less common you need the extra include.
>
> Best
>
> Bruno
>
> On Friday, May 3, 2019 at 4:07:56 PM UTC-4, Hamed Babaei wrote:
>>
>> Hello,
>>
>> I am planning to use a symmetric tensor rank-6 for third-order elastic 
>> constants in a nonlinear elasticity code. The tensor is supposed to be 
>> symmetric within each three pair indices as well as all possible orders of 
>> pair permutations. 
>> Namely, C_ijklmn = C_jiklmn = C_ijlkmn = C_ijklnm = C_klijmn = C_mnijkl 
>> ,.This tensor without any material symmetry will have 56 independent 
>> elements. 
>>
>> The SymmetricTensor class provides the symmetry within the pairs without 
>> the symmetry obtained by permutation of pairs and it will be totally fine 
>> and I can take care of the extra symmetry I need by writing the equal 
>> elements as many times as required.  
>>
>> However, the problem is that I can not even initialize a rank-6 tensor. 
>> Interestingly it works for rank-4 though. It seems that I am missing 
>> something in the instruction of the tensor. Only the following simplest 
>> line of code gives me an error.
>>
>> $ *SymmetricTensor<6, dim> C; *   
>>
>>
>> The error is:
>>
>> /home/hbabaei/deal.ii-candi/deal.II-v9.0.1/include/deal.II/base/symmetric_tensor.h:563:29:
>>  
>> error: incomplete type 
>> ‘dealii::internal::SymmetricTensorAccessors::StorageType<6, 3, double>’ 
>> used in nested name specifier
>>static const unsigned int n_independent_components
>>  ^
>> /home/hbabaei/deal.ii-candi/deal.II-v9.0.1/include/deal.II/base/symmetric_tensor.h:1046:1:
>>  
>> error: size of array is not an integral constant-expression
>>  SymmetricTensor::SymmetricTensor (const Number 
>> () [n_independent_components])
>>  ^
>> /home/hbabaei/deal.ii-candi/deal.II-v9.0.1/include/deal.II/base/symmetric_tensor.h:1046:100:
>>  
>> error: size of array is not an integral constant-expression
>>  SymmetricTensor::SymmetricTensor (const Number 
>> () [n_independent_components])
>>   
>>   ^
>> /home/hbabaei/deal.ii-candi/deal.II-v9.0.1/include/deal.II/base/symmetric_tensor.h:861:61:
>>  
>> error: invalid use of incomplete type ‘struct 
>> dealii::internal::SymmetricTensorAccessors::StorageType<6, 3, double>’
>>typedef typename base_tensor_descriptor::base_tensor_type 
>> base_tensor_type;
>>  ^
>> /home/hbabaei/deal.ii-candi/deal.II-v9.0.1/include/deal.II/base/symmetric_tensor.h:188:12:
>>  
>> error: declaration of ‘struct 
>> dealii::internal::SymmetricTensorAccessors::StorageType<6, 3, double>’
>>  struct StorageType;
>> ^
>>
>> It seems that I need to specify n_independent_components, but I do not 
>> know how and where I can do so.
>>
>> I appreciate it if you could give me any clue how to implement a rank-6 
>> tensor successfully.
>>
>> Thanks and regards,
>>
>

-- 
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] SymmetricTensor rank 6 implementation issue

[Hamed Babaei]

Hello,

I am planning to use a symmetric tensor rank-6 for third-order elastic 
constants in a nonlinear elasticity code. The tensor is supposed to be 
symmetric within each three pair indices as well as all possible orders of 
pair permutations. 
Namely, C_ijklmn = C_jiklmn = C_ijlkmn = C_ijklnm = C_klijmn = C_mnijkl 
,.This tensor without any material symmetry will have 56 independent 
elements. 

The SymmetricTensor class provides the symmetry within the pairs without 
the symmetry obtained by permutation of pairs and it will be totally fine 
and I can take care of the extra symmetry I need by writing the equal 
elements as many times as required.  

However, the problem is that I can not even initialize a rank-6 tensor. 
Interestingly it works for rank-4 though. It seems that I am missing 
something in the instruction of the tensor. Only the following simplest 
line of code gives me an error.

$ *SymmetricTensor<6, dim> C; *   


The error is:

/home/hbabaei/deal.ii-candi/deal.II-v9.0.1/include/deal.II/base/symmetric_tensor.h:563:29:
 
error: incomplete type 
‘dealii::internal::SymmetricTensorAccessors::StorageType<6, 3, double>’ 
used in nested name specifier
   static const unsigned int n_independent_components
 ^
/home/hbabaei/deal.ii-candi/deal.II-v9.0.1/include/deal.II/base/symmetric_tensor.h:1046:1:
 
error: size of array is not an integral constant-expression
 SymmetricTensor::SymmetricTensor (const Number () 
[n_independent_components])
 ^
/home/hbabaei/deal.ii-candi/deal.II-v9.0.1/include/deal.II/base/symmetric_tensor.h:1046:100:
 
error: size of array is not an integral constant-expression
 SymmetricTensor::SymmetricTensor (const Number () 
[n_independent_components])

^
/home/hbabaei/deal.ii-candi/deal.II-v9.0.1/include/deal.II/base/symmetric_tensor.h:861:61:
 
error: invalid use of incomplete type ‘struct 
dealii::internal::SymmetricTensorAccessors::StorageType<6, 3, double>’
   typedef typename base_tensor_descriptor::base_tensor_type 
base_tensor_type;
 ^
/home/hbabaei/deal.ii-candi/deal.II-v9.0.1/include/deal.II/base/symmetric_tensor.h:188:12:
 
error: declaration of ‘struct 
dealii::internal::SymmetricTensorAccessors::StorageType<6, 3, double>’
 struct StorageType;
^

It seems that I need to specify n_independent_components, but I do not know 
how and where I can do so.

I appreciate it if you could give me any clue how to implement a rank-6 
tensor successfully.

Thanks and regards,

-- 
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] Trouble using VectorTools::interpolate_based_on_material_id

[Hamed Babaei]

Dear Wolfgang,

In the meanwhile that the VectorTools::interpolate_based_on_material_id 
function is taking care of, I just did the same through the following code 
to change my solution vector in the specific regions with designated 
material_id, which may be used by those who have the same issue with this 
function until it is fixed. Hope it helps!
 
///
FEValues fe_values (fe, qf_cell, update_values  );
std::vector local_dof_indices(dofs_per_cell);

typename DoFHandler::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler_c.end();
for (; cell!=endc; ++cell)
  if (cell->is_locally_owned())
   {
  fe_values.reinit(cell);
  if (cell->material_id()==1)
   {
cell->get_dof_indices (local_dof_indices);
   for (unsigned int i=0; ihttp://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] Trouble using VectorTools::interpolate_based_on_material_id

[Hamed Babaei]

Dear Wolfgang: Thank you for your time addressing this issue.

-- 
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] Trouble using VectorTools::interpolate_based_on_material_id

[Hamed Babaei]

Dear Wolfgang,

Thank you for opening the Github issue. I've been looking around the 
function implementation and will be trying to figure out the issue. I will 
update you once I have some understanding.

Regards,

-- 
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] Trouble using VectorTools::interpolate_based_on_material_id

[Hamed Babaei]

I forgot to mention that I am using deal.ii 9.0.0 version. 
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.

Re: [deal.II] Trouble using VectorTools::interpolate_based_on_material_id

[Hamed Babaei]

Dear Wolfgang,
 

> In other words, there are no NaNs for me. What do you get? And what 
> version of 
> deal.II are you using? 
>
 
Thank you very much for your help. In fact, what you got is for the case 
that the material_id for the entire domain (both half of cubes) is included 
in function_map. My domain is a cube divided into two equal part at the 
plane x=0.5. One half has material_id=0 and the other half has 
material_id=1. So if you comment the following lines in Run:

const ConstantFunction constant_function_2(2.0); 

function_map[0] = _function_2;  

the half with material_id=0 has not included in the function_map and 
consequently, it leads to the -nan values.

P.S. for the dst vector elements, I considered constant value 1 and 2 for 
cells with material_id=1 and 0 respectively.

Thanks and regards, 

-- 
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: Trouble using VectorTools::interpolate_based_on_material_id

[Hamed Babaei]

I just realized I made a mistake and uploaded the executable instead of the 
mini code itself!  It is uploaded again here.

On Monday, December 10, 2018 at 11:17:40 AM UTC-6, Hamed Babaei wrote:
>
> Hello, 
>
> I would like to forcefully determine the solution at some regions of my 
> domain in which I solve an initial value problem. Therefore, I chose to 
> use VectorTools::interpolate_based_on_material_id function to implement the 
> desired constant value for solution vector within particular regions for 
> which I designate different material_ID. However, the problem is that 
> despite the function description saying: 
>
> "If a material id of some group of cells is missed in function_map, then 
> dst will not be updated in the respective degrees of freedom of the 
> output vector For example, if dst was successfully initialized to capture 
> the degrees of freedom of the dof_handler of the problem with all zeros 
> in it, then those zeros which correspond to the missed material ids will 
> still remain in dst even after calling this function."
>
> I receive -nan value for the cells whose material id is not included in 
> function_map when using the interpolate_based_on_material_id function.  
> However, if I include all the material ids in the function_map it works as 
> expected for the entire domain. This is not what I need because I would 
> like to keep solution constant only for a small part of the domain, not the 
> entire domain.
>
> I really appreciate it if you could give me any clue about where I am 
> making mistake using this function. Please find attached the minimal code 
> for your review. 
>
> Thanks and regards  
>
>

-- 
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.
#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 PhaseField
{
  using namespace dealii;

  typedef TrilinosWrappers::MPI::Vector vectorType;
  typedef TrilinosWrappers::SparseMatrix matrixType;


///

  template 
  class Solid
  {
  public:
Solid();

virtual
~Solid();

void
run();

  private:

voidmake_grid();

MPI_Comm mpi_communicator;
parallel::distributed::Triangulation triangulation;
ConditionalOStream   pcout;
FE_Qfe;
DoFHandler  dof_handler;
IndexSet locally_owned_dofs;
IndexSet locally_relevant_dofs;
  };


// constructor
  template 
  Solid::Solid()
:
mpi_communicator (MPI_COMM_WORLD),
triangulation (mpi_communicator,
   typename Triangulation::MeshSmoothing
   (Triangulation::smoothing_on_refinement |
Triangulation::smoothing_on_coarsening)),

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

//destructor
  template 
  Solid::~Solid()
  {
dof_handler.clear();
  }


///
  template 
  void Solid::make_grid()
  {

std::vector< unsigned int > repetitions(dim, 4);
if (dim == 3)
repetitions[dim-2] = 4;
repetitions[dim-3] = 4;

GridGenerator::subdivided_hyper_rectangle(triangulation,
	  repetitions,
  Point(0.0, 0.0, 0.0),
	  Point(1.0, 1.0, 1.0),
  	   true);


	for (typename Triangulation::active_cell_iterator
	cell = triangulation.begin_active();
	cell != triangulation.end(); ++cell)

		if(cell->center()[0]<0.5)
	  cell->set_material_id (1);
		else
		 cell->set_material_id (0);


  }

  //
  template 
  void Solid::run()
  {
make_grid();
dof_handler.distribute_dofs (fe);

locally_owned_dofs = dof_handler.locally_owned_dofs ();
vectorType  dst(locally_owned_dofs, mpi_communicator);

pcout << " Initial values: " << std::endl;
dst.print(std::cout);

std::map< types::material_id, const Function* > function_map;

 

[deal.II] Re: GridTools::collect_periodic_faces has problem with meshes generated by Gmsh or Abaqus when running in parallel.

[Hamed Babaei]

Dear Dr. Arndt,

Thank you very much for your help. The problem was resolved applying your 
comment. I was setting the indicators running over locally owned cells not 
all the cells.

Best regards,
Hamed 

On Saturday, February 17, 2018 at 11:20:32 AM UTC-6, Daniel Arndt wrote:
>
> Hamed,
>
> for collecting the periodic boundary pairs only the coarsest level is 
> important. 
> In particular, you have to make sure that you set the boundary indicators 
> correctly for all cells on this level.
> For the coarsest level every MPI process knows about all the cells since 
> the mesh is the same for all processes. 
> Summarizing: Simply run over all the cells when setting the boundary 
> indicators and not only the ones that are locally owned.
>
> 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] GridTools::collect_periodic_faces has problem with meshes generated by Gmsh or Abaqus when running in parallel.

[Hamed Babaei]

Hi All,

I'm having problem when applying periodicity on meshes generated with Gmsh 
or Abaqus. The problem arises when I run in parallel (in serial it works) 
getting the following error:

 void 
dealii::GridTools::match_periodic_face_pairs(std::set >&, std::set >&, int, 
std::vector&, const 
dealii::Tensor<1, typename FaceIterator::AccessorType:: space_dimension>&, 
const dealii::FullMatrix&) [with CellIterator = 
dealii::TriaIterator >; typename 
dealii::identity::type = 
dealii::TriaIterator >; typename 
FaceIterator::AccessorType = dealii::CellAccessor<3, 3>]
The violated condition was:
n_matches == pairs1.size() && pairs2.size() == 0
The name and call sequence of the exception was:
ExcMessage ("Unmatched faces on periodic boundaries")
Additional Information:
Unmatched faces on periodic boundaries 

I suspect that match_periodic_face_pairs function does not have access to 
the whole faces in parallel. 

I am wondering if you have encountered this issue before, or any clue what 
is the reason for this issue.

Thanks in advance.

Hamed

-- 
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: Multiple Triangulations

[Hamed Babaei]

Dear Dr. Bangert,

Since you have hanging nodes, do you apply hanging node constraints to 
> the displacement field you then put into MappingQEulerian? 
>

It turns out that I have not applied hanging node constraints properly. I 
just did so and the problem was solved.
Thank you very much for your help.

Best regards,
Hamed  

-- 
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: Adding periodicity to a thin sample with one element in thickness direction

[Hamed Babaei]

Hello all,

You are right. I was running in debug mode, slowing down the code.
Thank you very much for your incredible help. I appreciate your kindness 
and generosity to support the community. 

Best regards,
Hamed

-- 
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: Adding periodicity to a thin sample with one element in thickness direction

[Hamed Babaei]

Hello all,

I am so sorry, Wolfgang, if I failed to express my thoughts in a clear and 
proper manner. Let me restate my reply to J-P's comment, Hope it is better 
this time !
 

> I've attached a slightly amended version of your code that I think will 
> highlight what the perceived problem is. When you use GridGenerator 
> ::hyper_rectangle you start out with a single cell that is refined globally 
> 3 times, which in dim=3 results in 1*(2^3)^3 = 512 cells. For your other 
> grid you start off with 30*30*1 = 900 cells on the coarsest level, which 
> upon refinement becomes 900*(2^3)^3 = 460 800 cells. I think that you might 
> find that its taking a little while to produce the final refined grid.
>
 
when using subdivided_hyper_rectangle and starting with 30*30*1 = 900 cells 
I do not need to refine any more. Triangulation.refine_global () is 
required for the case of hyper_rectangle. 

However, I think it is not the producing of the final refined mesh that 
hangs the code but the add_periodicity when repetition parameter is 
relatively large(100 for instance). Since, it has been mentioned in the 
documentation that refinement of the course mesh should be applied after 
adding periodicity, I think the problem is that using repetition in 
subdivided_hyper_rectangle, I am refining the course mesh before 
add_periodicity.   
 

> Can you confirm that the code I've attached reaches the end, i.e. that you 
> see this result:
>
> $ mpirun -np 2 ./distib_tria_periodicity
>> Adding periodicity...  done.
>> Refinement iteration: 1
>> Refinement iteration: 2
>> Refinement iteration: 3
>> Finished.
>
>  
The code reaches the end for small repetition parameter and hangs for large 
ones.

I am wondering if I am right in thinking about the source of the problem. 

Thank you all for your time and prompt attention.  

Regards,
Hamed

-- 
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: Adding periodicity to a thin sample with one element in thickness direction

[Hamed Babaei]

Dear Jean-Paul,

Thank you for your help. For small number of repetition argument, for 
example less than 50, the code reaches the end. 
However, I just checked 100 for repetition and it has been for about an 
hour that the code is still stuck at the 
"triangulation.add_periodicity(periodicity_vector);" !!
It seems that repetition factor refines the mesh and makes add_periodicity 
too slow.

Thanks


On Tuesday, May 16, 2017 at 12:03:14 PM UTC-5, Jean-Paul Pelteret wrote:
>
> Dear Hamed,
>
> I could not find a bug in the code you provided. It worked perfectly fine 
> for me using deal.II 8.5 and with the number of MPI processes = {1,2}. 
>
> I've attached a slightly amended version of your code that I think will 
> highlight what the perceived problem is. When you use GridGenerator 
> ::hyper_rectangle you start out with a single cell that is refined globally 
> 3 times, which in dim=3 results in 1*(2^3)^3 = 512 cells. For your other 
> grid you start off with 30*30*1 = 900 cells on the coarsest level, which 
> upon refinement becomes 900*(2^3)^3 = 460 800 cells. I think that you might 
> find that its taking a little while to produce the final refined grid.
>
> Can you confirm that the code I've attached reaches the end, i.e. that you 
> see this result:
>
> $ mpirun -np 2 ./distib_tria_periodicity
>> Adding periodicity...  done.
>> Refinement iteration: 1
>> Refinement iteration: 2
>> Refinement iteration: 3
>> Finished.
>
>
> Regards,
> Jean-Paul
>
> On Tuesday, May 16, 2017 at 6:17:07 PM UTC+2, Hamed Babaei wrote:
>>
>> Dear Jean-Paul and Daniel,
>>
>> I provided a minimal code which demonstrates the issue. There is no 
>> problem when I use  GridGenerator::hyper_rectangle to create a course mesh 
>> and after adding periodicity, triangulation.refine_global is called.
>> However, when I use subdivided_hyper_rectangle with repetition argument, 
>> code gets stuck ! which I quess is because of repetition argument that 
>> refines the mesh before adding periodicity. 
>>
>> //
>> #include 
>> #include 
>> #include 
>> #include 
>> #include 
>> #include 
>> #include 
>>
>> #include 
>> #include 
>>
>> #include 
>> #include 
>> /
>> using namespace dealii;
>>
>> template 
>>   class Solid
>>   {
>>   public:
>> Solid();
>>
>> voidmake_grid();
>>
>>   private:
>> MPI_Comm mpi_communicator;
>> parallel::distributed::Triangulation triangulation;
>>
>>   };
>>
>> template 
>>   Solid::Solid()
>> :
>> mpi_communicator (MPI_COMM_WORLD),
>> triangulation (mpi_communicator,
>>typename Triangulation::MeshSmoothing
>>(Triangulation::smoothing_on_refinement |
>> Triangulation::smoothing_on_coarsening))
>>
>> {}
>>
>> template 
>>   void Solid::make_grid()
>> {
>>   std::vector< unsigned int > repetitions(dim, 30);
>>   if (dim == 3)
>>   repetitions[dim-1] = 1;
>>
>> //  GridGenerator::subdivided_hyper_rectangle(triangulation,
>> //repetitions,
>> //Point(0.0, 0.0, -0.5),
>> //Point(20.0, 20.0, 0.5),
>> //true);
>>
>>
>>GridGenerator::hyper_rectangle(triangulation,
>>   Point(0.0, 0.0, 0.0),
>>   Point(20.0, 20.0, 20.0),
>>   true);
>>
>>std::vector<GridTools::PeriodicFacePair> parallel::distributed::Triangulation::cell_iterator> >
>>periodicity_vector;
>>
>>GridTools::collect_periodic_faces(triangulation,
>>  /*b_id1*/ 0,
>>  /*b_id2*/ 1,
>>  /*direction*/ 0,
>>   periodicity_vector);
>>
>>   triangulation.add_periodicity(periodicity_vector);
>>
>>
>>   triangulation.refine_global (3);
>>
>> }
>>
>> int main (int argc, char *argv[])
>>   {
>> Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
>> Solid<3> solid_3d;
>> solid_3d.make_grid();
>>   }
>> //
>>
>> Any clue would be appreciated.
>>
>> Thanks and regards,
>> Ha

[deal.II] Re: Adding periodicity to a thin sample with one element in thickness direction

[Hamed Babaei]

Dear Jean-Paul and Daniel,

I provided a minimal code which demonstrates the issue. There is no problem 
when I use  GridGenerator::hyper_rectangle to create a course mesh and 
after adding periodicity, triangulation.refine_global is called.
However, when I use subdivided_hyper_rectangle with repetition argument, 
code gets stuck ! which I quess is because of repetition argument that 
refines the mesh before adding periodicity. 

//
#include 
#include 
#include 
#include 
#include 
#include 
#include 

#include 
#include 

#include 
#include 
/
using namespace dealii;

template 
  class Solid
  {
  public:
Solid();

voidmake_grid();

  private:
MPI_Comm mpi_communicator;
parallel::distributed::Triangulation triangulation;

  };

template 
  Solid::Solid()
:
mpi_communicator (MPI_COMM_WORLD),
triangulation (mpi_communicator,
   typename Triangulation::MeshSmoothing
   (Triangulation::smoothing_on_refinement |
Triangulation::smoothing_on_coarsening))

{}

template 
  void Solid::make_grid()
{
  std::vector< unsigned int > repetitions(dim, 30);
  if (dim == 3)
  repetitions[dim-1] = 1;

//  GridGenerator::subdivided_hyper_rectangle(triangulation,
//repetitions,
//Point(0.0, 0.0, -0.5),
//Point(20.0, 20.0, 0.5),
//true);


   GridGenerator::hyper_rectangle(triangulation,
  Point(0.0, 0.0, 0.0),
  Point(20.0, 20.0, 20.0),
  true);

   std::vector<GridTools::PeriodicFacePair::cell_iterator> >
   periodicity_vector;

   GridTools::collect_periodic_faces(triangulation,
 /*b_id1*/ 0,
 /*b_id2*/ 1,
 /*direction*/ 0,
  periodicity_vector);

  triangulation.add_periodicity(periodicity_vector);


  triangulation.refine_global (3);

}

int main (int argc, char *argv[])
  {
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
Solid<3> solid_3d;
solid_3d.make_grid();
  }
//

Any clue would be appreciated.

Thanks and regards,
Hamed

On Monday, May 15, 2017 at 8:44:24 AM UTC-5, Jean-Paul Pelteret wrote:
>
> Dear Hamad,
>
> Would you mind please posting a minimal example that demonstrates the 
> issue? Since both Daniel and I are surprised that this is the behaviour 
> that you're experiencing, I would say that this is an unexpected outcome.
>
> Regards,
> Jean-Paul
>
> On Monday, May 15, 2017 at 3:40:20 PM UTC+2, Hamed Babaei wrote:
>>
>> Dear Daniel,
>>
>> As far as I see, GridGenerator::subdivided_hyper_rectangle only calls 
>>> "Triangulation::create_triangulation" and hence creates a coarse 
>>> triangulation.
>>> Can you confirm either?
>>>
>> If  GridGenerator::subdivided_hyper_rectangle creates a coarse 
>> triangulation, a question arises that when it uses repetition argument 
>> which specifies number of elements in every direction. 
>> Besides, when I use GridGenerator::subdivided_hyper_rectangle with 
>> repetition argument, the code gets stuck at add_periodicity and never ends. 
>>
>> Thanks and regards,
>> Hamed  
>>
>

-- 
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: Adding periodicity to a thin sample with one element in thickness direction

[Hamed Babaei]

Dear Daniel,

As far as I see, GridGenerator::subdivided_hyper_rectangle only calls 
> "Triangulation::create_triangulation" and hence creates a coarse 
> triangulation.
> Can you confirm either?
>
If  GridGenerator::subdivided_hyper_rectangle creates a coarse 
triangulation, a question arises that when it uses repetition argument 
which specifies number of elements in every direction. 
Besides, when I use GridGenerator::subdivided_hyper_rectangle with 
repetition argument, the code gets stuck at add_periodicity and never ends. 

Thanks and regards,
Hamed  

-- 
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: Adding periodicity to a thin sample with one element in thickness direction

[Hamed Babaei]

Dear Jean-Paul,

Thank you for your prompt help.

Although I don't really understand your specific question, I am guessing 
> that you were wondering if its possible to use the "repetitions" parameter 
> to bias the initial mesh gradation before adding periodicity (since the 
> periodic face pairs must be chosen on the on the coarsest mesh).
>
 
Yes, you are right. In fact, I do not want to refine the mesh but generate 
the initial mesh in away that mesh discritisation is done after adding 
periodicity.
 

> To the best of my knowledge, GridGenerator::subdivided_hyper_rectangle 
> generates a coarsest level mesh (i.e. all cells have level = 0) with the 
> given discretisation in each coordinate direction. 
>

It seems that GridGenerator::subdivided_hyper_rectangle needs repetition 
argument and complete mesh generation in one step. so cells have level=1 .
 

> You could also consider using anisotropic refinement (see step-30 
> ), but this 
> feature is not so commonly used (as far as I can tell) and it is not 
> supported by all Triangulations 
> 
> .
>

Since my code is parallel, I am wondering what is the alternative for me.


Thanks and regards,
Hamed
 

-- 
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] Adding periodicity to a thin sample with one element in thickness direction

[Hamed Babaei]

Hello,

I want to apply periodic boundary condition on a 3D thin square plate with 
one element in the thickness direction. I have already used the following 
for meshing

std::vector< unsigned int > repetitions(dim, 100);
if (dim == 3)
repetitions[dim-1] = 1;
GridGenerator::subdivided_hyper_rectangle(triangulation,
   repetitions,
  
 /*bottom_left*/ Point(0.0, 0.0, -0.5),
  /* 
top_right*/ Point(20.0, 20.0, 0.5),
  true); 

However, when applying periodic condition, mesh refinement should be 
applied after doing "triangulation.add_periodicity(periodicity_vector)".
So I have to first just use "hyper_rectangle" without any repetition 
argument to create the geometry, then after I added periodicity to the 
faces, mesh refinement can be applied.
Yet, I am wondering which function can apply different refinement for 
different directions, since " triangulation.refine_global ()" do the same 
for all directions.

Thanks and regards,
Hamed

-- 
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 output a single scalar in a parallel code

[Hamed Babaei]

Dear Rajat,
 

> ... which I then write in a file on a master process.
>

I was wondering how to write in a file on a master process so that there 
would be just one output file not as many as processor exist. I use the 
following commands to write in my file in serial code :

std::ofstream myfile;
myfile.open ("resultant_strass.txt");
myfile<< resultant_stress<

[deal.II] Re: Periodic Boundary Condition

[Hamed Babaei]

 Dear Daniel,

You can obtain local DoF values from a global FE vector using 
> FEValuesBase::get_function_values [1], create the right-hand side locally 
> as you wish and use ConstraintMatrix::distribute_local_to_global to create 
> the global right-hand side. This is outlined in step-21 for example. 
>

Considering for example the term  [*M*+k*theta**A*]**U* in my right-hand 
side, since the mass matrix(M) and laplace matrix(A) in the local level 
have the size of *dofs_per_cell * dofs_per_cell* I need to have a local 
solution vector to multiply with them. However, 
FEValuesBase::get_function_values returns solution values at quadrature 
points in a vector with the size of  *n_q_points*. I was wondering if there 
is a function that returns nodal solution values for every cell.

Thank you very much for your time.

Best regards,
Hamed

-- 
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: Periodic Boundary Condition

[Hamed Babaei]

Dear Daniel,

Happy new year. I wish a wonderful new year for you !

This is only (directly) possible if you are not using inhomogeneous 
> constraints. In that case, you can call 
> ConstraintMatrix::distribute_local_to_global separately for matrix and 
> right-hand side. 
>

I am struggling to apply homogeneous periodic BC constraints in assembly 
since ConstraintMatrix::condense is not usable in parallel and I can not 
use ConstraintMatrix::distribute_local_to_global since I construct 
system_matrix and system_rhs in global level as a linear combination of 
pre-costructed mass_matrix(M), laplace_matrix(A), nonlinear 
term(nl_term),nonlinear matrix(nl_matrix), solution (U) and 
old_solution(U_old) vectors as follows:

system_rhs= [M+k*theta*A]*U-[M-k*(1-theta)*A]*U_old -k*nl_term(U,U_old) 
 in which theta is Crank_Nicolson parameter 
system_matrix= M+k*theta*A-k*theta*nl_matrix(U, U_old)

However, to be able to use ConstraintMatrix::distribute_local_to_global, I 
can assemble system_matrix in local level instead of global level. Having 
that said, I can not do so for the system_rhs because solution and 
old_solution vectors are global vectors. 

Any hint regarding how to apply constraint on my system_rhs would be 
appreciated.

Thanks,
Hamed



-- 
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] Alternative for ConstraintMatrix::condense in parallel codes?

[Hamed Babaei]

Hi all,

I have used* constraints.condense* function in a serial code to apply 
periodic constraints on system_matrix and system_rhs after their assembly 
in which I do not use *constraints.distribute_local_to_global* for 
cell_matrix and cell_rhs but I first assemble mass_matrix , laplace_matrix, 
nl_matrix and nl_term (the same approach as in step-25) and then 
syetem_matrix and system_rhs are obtained after some manipulations on the 
mentioned assembled matrices and vectors.

The problem is that* constraints.condense* doesn't work with 
*TrilinosWrappers::SparseMatrix 
*so the question is when we can not use *constraints.distribute_local_to_global 
*how we can apply constraints in the assembly for parallel codes. In the 
https://www.dealii.org/8.4.1/doxygen/deal.II/group__constraints.html 

 the 
following is mentioned which I can't truly understand. 

"The condensation functions exist for different argument types: 
SparsityPattern 
, 
SparseMatrix 
 and 
BlockSparseMatrix 
. 
Note that there are no versions for arguments of type 
PETScWrappers::SparseMatrix() 

 or 
any of the other PETSc or Trilinos matrix wrapper classes. This is due to 
the fact that it is relatively hard to get a representation of the sparsity 
structure of PETSc matrices, and to modify them efficiently; this holds in 
particular, if the matrix is actually distributed across a cluster of 
computers. If you want to use PETSc/Trilinos matrices, you can either copy 
an already condensed deal.II matrix, or assemble the PETSc/Trilinos matrix 
in the already condensed form, see the discussion below."

Thanks and regards,
Hamed

-- 
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: Periodic boundary conditions seem to be not applied

[Hamed Babaei]

Dear Daniel,

I just resolved the problem. Because constraints.close function resolves 
chains of constraints, It interfered with set_inhomogeinity process. So I 
had to issue constraints.close before applying inhomogeinities.

Thank you very much for your helpful hints. I appreciate it.

Hamed  

-- 
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: Periodic boundary conditions seem to be not applied

[Hamed Babaei]

Dear Daniel,
 

> It seems that you are now constraining all the components, i.e. 
> u_x(0,y)=u_x(L,y) and u_y(0,y)=u_y(L,y).
> Is this really what you want to do? Otherwise, you should have a 
> ComponentMask in your call to 
> DoFTools::make_periodicity_constraints as well, i.e.
>
>  You are right Daniel. I want to constraint all displacement components. 
If I constraint U_x for left and right faces and U_y for top and bottom 
faces, I do not confront the problem I have now. The problem is that the 
DoFs residing at the bottom-left node (where two faces on which we apply 
inhomogenity, interface) the inhomogenity for the U_y is doubled despite 
the similarity of the code for x and y direction !!!  However, this happens 
by "constraints.close" function since when I uncomment it, the constraints 
are exactly what is coded. Nevertheless, the constraints.close is required 
and can not be omitted.  

Again, I don't see an immediate reason why x- and y-direction don't behave 
> similarly. Can you provide a minimal working code showing this?
>

The minimal code is attached. It would be appreciated if you could 
investigate what is wrong!

Thanks,
Hamed 
 

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


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

#include 
#include 

namespace PhaseField
{
  using namespace dealii;

  template 
  class Solid
  {
  public:
Solid();

virtual
~Solid();

void
run();

private:

voidmake_grid();
voidmake_constraints();


Triangulation   triangulation;
const FESystem  fe;
DoFHandler  dof_handler;
ConstraintMatrix constraints;

  };

  template 
  Solid::Solid()
:

fe(FE_Q(1), dim),
dof_handler(triangulation)
  {}

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

  template 
  void Solid::make_grid()
  {
GridGenerator::hyper_cube (triangulation, -1, 1,true);
dof_handler.distribute_dofs(fe);
  }

 
  template 
void Solid::make_constraints()
{
 constraints.clear();

  const FEValuesExtractors::Scalar x_displacement(0);
  const FEValuesExtractors::Scalar y_displacement(1);
  const FEValuesExtractors::Scalar z_displacement(2);

   DoFTools::make_periodicity_constraints(dof_handler,
  	  /*b_id*/ 0,
  	  /*b_id*/ 1,
  	  /*direction*/ 0,
  	  constraints);

   DoFTools::make_periodicity_constraints(dof_handler,
  	  /*b_id*/ 2,
  	  /*b_id*/ 3,
  	  /*direction*/ 1,
  	  constraints);





{
  IndexSet selected_dofs_x;
  std::set< types::boundary_id > boundary_ids_x= std::set();
  boundary_ids_x.insert(0);

  DoFTools::extract_boundary_dofs(dof_handler,
 fe.component_mask(x_displacement),
			 selected_dofs_x,
			 boundary_ids_x);
  unsigned int nb_dofs_face_x = selected_dofs_x.n_elements();
  IndexSet::ElementIterator dofs_x = selected_dofs_x.begin();

 for(unsigned int i = 0; i < nb_dofs_face_x; i++)
  {
	constraints.set_inhomogeneity(*dofs_x, 0.1);
  dofs_x++;
  }
}

{
IndexSet selected_dofs_y;
std::set< types::boundary_id > boundary_ids_y= std::set();
boundary_ids_y.insert(2);

DoFTools::extract_boundary_dofs(dof_handler,
   fe.component_mask(y_displacement),
  			   selected_dofs_y,
  			   boundary_ids_y);
unsigned int nb_dofs_face_y = selected_dofs_y.n_elements();
IndexSet::ElementIterator dofs_y = selected_dofs_y.begin();

for(unsigned int i = 0; i < nb_dofs_face_y; i++)
{

  	constraints.set_inhomogeneity(*dofs_y, 0.3 );
dofs_y++;
}
  }




  constraints.close();

}

  //
 template 
 void Solid::run()
  {
   make_grid();
   make_constraints();
   constraints.print(std::cout);

  }

}

[deal.II] Re: Periodic boundary conditions seem to be not applied

[Hamed Babaei]

Dear Daniel and Timo,

Thank you very much for your advises. I have figured out  the problem with 
the periodic constraints but I do not know where this problem comes from.
Let me explain the issue more. As I mentioned earlier, I have applied 
periodic BC for displacement field and I want to apply boundary 
displacement as well using constraints.set_inhomogeneit. I mean, 
u(0,y)=u(L,y)+delta_x  for the left and right faces and 
u(x,0)=u(x,L)+delta_y for the bottom and top faces. 

 The following is my periodic constraints and the code for boundary 
displacement. 
//
template 
void Solid::make_constraints(const int _nr)
{
  if (it_nr > 1)
return;

  constraints.clear();

  const bool apply_dirichlet_bc = (it_nr == 0);
  

   DoFTools::make_periodicity_constraints(dof_handler,
/*b_id*/ 0,
/*b_id*/ 1,
/*direction*/ 0,
constraints);
   DoFTools::make_periodicity_constraints(dof_handler,
/*b_id*/ 2,
/*b_id*/ 3,
/*direction*/ 1,
constraints);

// This block add to the periodicity constraint the boundary 
displacement
   // in x direction for left and right faces and in y direction 
for top and bottom faces


{
{
  IndexSet selected_dofs_x;
  std::set< types::boundary_id > boundary_ids_x= 
std::set();
  boundary_ids_x.insert(0);

  DoFTools::extract_boundary_dofs(dof_handler,

fe.component_mask(x_displacement),
selected_dofs_x,
boundary_ids_x);
  unsigned int nb_dofs_face_x = selected_dofs_x.n_elements();
  IndexSet::ElementIterator dofs_x = selected_dofs_x.begin();

 
  for(unsigned int i = 0; i < nb_dofs_face_x; i++)
  {
 constraints.set_inhomogeneity(*dofs_x, (apply_dirichlet_bc ? 
5e-2 : 0.0));
  dofs_x++;
  }
}

{
IndexSet selected_dofs_y;
std::set< types::boundary_id > boundary_ids_y= 
std::set();
boundary_ids_y.insert(2);

DoFTools::extract_boundary_dofs(dof_handler,

 fe.component_mask(y_displacement),
   selected_dofs_y,
   boundary_ids_y);
unsigned int nb_dofs_face_y = selected_dofs_y.n_elements();
IndexSet::ElementIterator dofs_y = selected_dofs_y.begin();

for(unsigned int i = 0; i < nb_dofs_face_y; i++)
{

  constraints.set_inhomogeneity(*dofs_y, (apply_dirichlet_bc ? 
-5e-2: 0.0));
dofs_y++;
}
  }
}


//and here I fix one displacement DoF in every direction  
  }
  }

  constraints.close();

}

//
I printed the constraint matrix for a domain of just one 2D linear element 
(4 node and 8 DoFs) as follows:  
   
   4,5___6,7
 |  | 
 |  |
 |  | 
 |  | 
   0,1__ 2,3 

when I just apply delta_x:

0 6: 1

0: 0.05

1 7: 1

2 6: 1

3 7: 1

4 6: 1

4: 0.05

5 7: 1


when I just have delta_y:

0 6: 1

1 7: 1

1: -0.1

2 6: 1

3 7: 1

3: -0.05

4 6: 1

5 7: 1

As you can see, the inhomogenity applied on DoF 1 is ironically doubled and 
that's why the corner element get distorted.

I was wondering why despite the identical code 
for constraints.set_inhomogeneity for the x and y direction, I have this 
problem.

Any hint would be appreciated,

Thanks,
Hamed

-- 
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: Periodic boundary conditions seem to be not applied

[Hamed Babaei]

Dear Daniel,
 

> I would suggest you to look at the constraints in the ConstraintMatrix 
> using ConstraintMatrix::print[1] for a mesh consisting of has few cells as 
> possible.
> Using DoFTools::map_dofs_to_support_points[2] you can then verify that 
> these are correct.
>

Issuing " constraints.print(std::cout);" I printed my Constraint matrix but 
I have two problems regarding using 
DoFTools::map_dofs_to_support_points(const Mapping< dim, spacedim > 
, const DoFHandler< dim, spacedim > _handler, std::vector< 
Point< spacedim > > _points); . First, I do not know how to 
construct and initialize the mapping argument since I could not find the 
constructor in the documentation. Second, I do not know how to visualize 
the result of this object. 

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.

[deal.II] Re: Periodic boundary conditions seem to be not applied

[Hamed Babaei]

I just realized that the problem comes from that part of the code that 
applies relative displacement between bottom and top faces, since there is 
not corner element distortion when having u_x (0,y)=u_x(L,y)+delta_x and I 
uncomment the second part of the code namely,

 {
IndexSet selected_dofs_y;
std::set< types::boundary_id > boundary_ids_y= 
std::set();
boundary_ids_y.insert(2);

DoFTools::extract_boundary_dofs(dof_handler,

 fe.component_mask(y_displacement),
selected_dofs_y,
 boundary_ids_y);
unsigned int nb_dofs_face_y = selected_dofs_y.n_elements();
IndexSet::ElementIterator dofs_y = selected_dofs_y.begin();

double relative_displacement_y;
if (timestep <10)
  relative_displacement_y = 0.05e-9;
else
 relative_displacement_y = 0.005e-9;

for(unsigned int i = 0; i < nb_dofs_face_y; i++)
{
   constraints.set_inhomogeneity(*dofs_y, (apply_dirichlet_bc ? 
relative_displacement_y : 0.0));
dofs_y++;
}

I cant realize what is the difference between relative displacement 
implementation in x and y direction that the code works for x direction but 
not for y direction while for both the code is identical !!!

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.

[deal.II] Re: Periodic boundary conditions seem to be not applied

[Hamed Babaei]

Dear Daniel,

I have applied periodic boundary condition on displacement field for a 
square 2D domain and in order to apply boundary displacement I want to 
consider relative displacement between pair faces, namely, 
u_x (0,y)=u_x(L,y)+delta_x  and u_y (x,0)=u_y(x,L)+delta_y , using 
"constraints.set_inhomogeneity" 
. I have used the following code for relative displacement between pair 
faces:

 {
  IndexSet selected_dofs_x;
  std::set< types::boundary_id > boundary_ids_x= 
std::set();
  boundary_ids_x.insert(0);

  DoFTools::extract_boundary_dofs(dof_handler,
  
 fe.component_mask(x_displacement),
   selected_dofs_x,
   boundary_ids_x);

  unsigned int nb_dofs_face_x = selected_dofs_x.n_elements();
  IndexSet::ElementIterator dofs_x = selected_dofs_x.begin();

  double relative_displacement_x;
  if (timestep <10)
 relative_displacement_x = -0.1e-9;
  else
 relative_displacement_x = -0.01e-9;

  for(unsigned int i = 0; i < nb_dofs_face_x; i++)
  {
constraints.set_inhomogeneity(*dofs_x, (apply_dirichlet_bc ? 
relative_displacement_x : 0.0));
  dofs_x++;
  }
}

{
IndexSet selected_dofs_y;
std::set< types::boundary_id > boundary_ids_y= 
std::set();
boundary_ids_y.insert(2);

DoFTools::extract_boundary_dofs(dof_handler,

 fe.component_mask(y_displacement),
   selected_dofs_y,
boundary_ids_y);
unsigned int nb_dofs_face_y = selected_dofs_y.n_elements();
IndexSet::ElementIterator dofs_y = selected_dofs_y.begin();

double relative_displacement_y;
if (timestep <10)
  relative_displacement_y = 0.05e-9;
else
   relative_displacement_y = 0.005e-9;

for(unsigned int i = 0; i < nb_dofs_face_y; i++)
{
  constraints.set_inhomogeneity(*dofs_y, (apply_dirichlet_bc ? 
relative_displacement_y : 0.0));
dofs_y++;
}
  }

Periodic boundary condition works great without boundary displacement and I 
can apply boundary displacement without any problem in two conditions : 
First, considering displacement between only two faces for example left and 
right faces, second, applying periodicity for only one components of 
displacement on pair faces, namely, u_x for left and right faces and u_y 
for top and bottom faces. However, considering periodicity for the both 
displacement components on pair faces and at the same time, applying both 
u_x (0,y)=u_x(L,y)+delta_x  and u_y (x,0)=u_y(x,L)+delta_y, the 
corner elements get distorted. In the attached picture you can see what 
happens.

Any hint about the source of error would be appreciated.

Thanks,
Hamed


-- 
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 apply non-homogeneous Neumann boundary conditions when we do not have exact solution

[Hamed Babaei]

Hello all,

Following step-7 I am going to apply non-homogeneous Neumann boundary 
conditions. The problem is that, in step-7, Neumann value is computed using 
exact solution, however in real-life problems we do not know exact solution 
beforehand. 

Neumann value is computed in step-7 as follows:

const double neumann_value
= (exact_solution.gradient (fe_face_values.quadrature_point 
(q_point))
 
*
fe_face_values.normal_vector 

(q_point));

what I have added to my assembly to apply non-homogeneous Neumann boundary 
conditions follows in which I compute Newman value by red lines:

for (unsigned int face_number=0; 
face_numberface(face_number)->at_boundary()) { fe_face_values.reinit (cell, 
face_number); fe_face_values.get_function_gradients(solution, 
solution_gradients); for (unsigned int q_point=0; q_point of file 

 
in function
void 
dealii::SparseMatrix::add(dealii::SparseMatrix::size_type, 
dealii::SparseMatrix::size_type, number) [with number = double; 
dealii::SparseMatrix::size_type = unsigned int]
The violated condition was: 
dealii::numbers::is_finite(value)
The name and call sequence of the exception was:
ExcNumberNotFinite(std::complex(value))
Additional Information: 
In a significant number of places, deal.II checks that some intermediate 
value is a finite number (as opposed to plus or minus infinity, or NaN/Not 
a Number). In the current function, we encountered a number that is not 
finite (its value is (nan,0) and therefore violates the current assertion.

This may be due to the fact that some operation in this function created 
such a value, or because one of the arguments you passed to the function 
already had this value from some previous operation. In the latter case, 
this function only triggered the error but may not actually be responsible 
for the computation of the number that is not finite.

There are two common cases where this situation happens. First, your code 
(or something in deal.II) divides by zero in a place where this should not 
happen. Or, you are trying to solve a linear system with an unsuitable 
solver (such as an indefinite or non-symmetric linear system using a 
Conjugate Gradient solver); such attempts oftentimes yield an operation 
somewhere that tries to divide by zero or take the square root of a 
negative value.

In any case, when trying to find the source of the error, recall that the 
location where you are getting this error is simply the first place in the 
program where there is a check that a number (e.g., an element of a 
solution vector) is in fact finite, but that the actual error that computed 
the number may have happened far earlier. To find this location, you may 
want to add checks for finiteness in places of your program visited before 
the place where this error is produced.One way to check for finiteness is 
to use the 'AssertIsFinite' macro.

Stacktrace:
---
#0  /home/hbabaei/deal.ii-candi/deal.II-v8.4.1/lib/libdeal_II.g.so.8.4.1: 
dealii::SparseMatrix::add(unsigned int, unsigned int, double)
#1  /home/hbabaei/deal.ii-candi/deal.II-v8.4.1/lib/libdeal_II.g.so.8.4.1: 
void 
dealii::ConstraintMatrix::distribute_local_to_global(dealii::FullMatrix const&, 
dealii::Vector const&, 
std::vector > const&, 
dealii::SparseMatrix&, dealii::Vector&, bool, 
dealii::internal::bool2type) const
#2  ./periodic-serial: PhaseField::Solid<3>::assemble_system_eta()
#3  ./periodic-serial: PhaseField::Solid<3>::solve_nonlinear_timestep_eta()
#4  ./periodic-serial: PhaseField::Solid<3>::run()
#5  ./periodic-serial: main




It would be appreciated if you could give any hint where I may be making 
mistake.

Thanks
Hamed

-- 
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: what to do with Neumann term in week formulation when applying Periodic Boundary Condition

[Hamed Babaei]

Another Point I would like to mention is that I checked constant values for 
neumann_value( for example, 0.7) and the solver converges. So I think 
problem comes from the way I compute neumann_value from old_solution.

I has followed step-7, in which non-homogeneous Neumann boundary values are 
applied using exact_solution for obtaining neumann_value. somewhere in 
documentation it is stated that:

"Finally we define an object denoting the exact solution function. We will 
use it to compute the Neumann values at the boundary from it. Usually, one 
would of course do so using a separate object, in particular since the 
exact solution is generally unknown while the Neumann values are 
prescribed. We will, however, be a little bit lazy and use what we already 
have in information. Real-life programs would to go other ways here, of 
course."

So the question is how to use FEM solution to obtain Neumann_values?

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.

[deal.II] what to do with Neumann term in week formulation when applying Periodic Boundary Condition

[Hamed Babaei]

Hi all,

I am considering Periodic Boundary Condition for a problem which is very 
simillar to step-25. Before applying periodic condition, I considered 
homogeneous Newmann Boundary condition in external surfaces by simply 
omitting Neumann term in week formulation (v , n.grad(u) )=0. However, 
since Neumann BC has meaning when we have external surface and now I have 
periodicity condition, It seems I can not have both periodic and neumann BC 
simultaneously. Therefore I am going to eliminate Neumann condition by 
keeping the term (v, n.grad(u)) in my right hand side as a generalized 
force term.  I added the following in my assembly to do so

 
  for (unsigned int face_number=0; 
face_numberface(face_number)->at_boundary())

   {
 fe_face_values_eta.reinit (cell, face_number);

 
fe_face_values_eta.get_function_gradients(old_solution_eta, 
solution_eta_gradients);

 for (unsigned int q_point=0; q_point

[deal.II] Re: Periodic boundary conditions seem to be not applied

[Hamed Babaei]

 Daniel,

You are completely right. The problem was resolved applying your comment.

Thank you very much,

Hamed

-- 
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: Periodic boundary conditions seem to be not applied

[Hamed Babaei]

Dear Bastian and Daniel,
 

> That's what I was assuming. So yesterday I did it with this code :
>   // This block add to the periodicity constraint the little 
> compression we want
>   {
> IndexSet selected_dofs_left;
> std::set< types::boundary_id > boundary_ids_left= 
> std::set();
> boundary_ids_left.insert(1);
> DoFTools::extract_boundary_dofs(dof_handler_ref,
>   fe.component_mask(x_displacement), selected_dofs_left, 
> boundary_ids_left);
> unsigned int nb_dofs_face_left = selected_dofs_left.n_elements();
> IndexSet::ElementIterator dofs_left = selected_dofs_left.begin();
> for(unsigned int i = 0; i < nb_dofs_face_left; i++)
> {
> //constraints.add_line(*dofs_left );
> constraints.set_inhomogeneity(*dofs_left, apply_dirichlet_bc ? 
> displacement_side_1 : 0.0);
> dofs_left++;
> }
>   }
> And it works well. Thanks a lot for your help.I think this topic can be 
> set to *Answered*!
>


I am going to do the same as Bastian, namely periodic boundary condition 
for displacement  such that u( 0, y ) = u( 1, y ) + *lambda ,*
* so I used the recommended above code by Bastian for applying inhomogenity 
to predefined periodicity, in order to have reletive dispacement between 
right and left faces,*
*but my results shows that both left and right faces are moving equally to 
the same direction, not relative to each other.*

*In addition, I activated the uncommented line, " *
 //constraints.add_line(*dofs_left );" otherwise I got this error: "call 
add_line() before calling set_inhomogeneity()"

This is my periodic constraints( left side boundry_id = 0  and right 
side boundry_id = 1)  :

 {
 std::vector >
 periodicity_vector_x;

GridTools::collect_periodic_faces(dof_handler,
/*b_id1*/ 0,
/*b_id2*/ 1,
/*direction*/ 0,
  periodicity_vector_x);
 DoFTools::make_periodicity_constraints 
(periodicity_vector_x,
   constraints,
   fe.component_mask(x_displacement));
  }

  {
  IndexSet selected_dofs_right;
  std::set< types::boundary_id > boundary_ids_right= 
std::set();
  boundary_ids_right.insert(1);

  DoFTools::extract_boundary_dofs(dof_handler,
  
 fe.component_mask(x_displacement),
selected_dofs_right,
boundary_ids_right);
  unsigned int nb_dofs_face_right = 
selected_dofs_right.n_elements();
  IndexSet::ElementIterator dofs_right = 
selected_dofs_right.begin();
  for(unsigned int i = 0; i < nb_dofs_face_right; i++)
  {
// constraints.add_line(*dofs_right );
 constraints.set_inhomogeneity(*dofs_right, (apply_dirichlet_bc 
? 0.01e-9 : 0.0));
  dofs_right++;
  }
  }

Any hint would be appreciated.

Thanks,
Hamed

-- 
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: Periodic Boundary Condition

[Hamed Babaei]

Dear Daniel,

As you can see my assembly in the first post, my system_matrix and 
system_rhs are not made directly from cell_system_matrix or 
cell_system_rhs, but they are the result of some manipulations on 
mass_matrix, laplase_matrix and nl_matrix (the same approach as step_25) 
that have already been assembled (in the red part of the code). 

To simplify the question, lets say we want to apply any sort of constraints 
in the assembly part of step-25. I was wondering how to do use 
ConstraintMatrix::distribute_local_to_global, while there is no 
cell_system_matrix and cell_system_rhs.

Thanks,
Hamed

-- 
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: Periodic Boundary Condition

[Hamed Babaei]

In addition I just tried to apply periodic boundary condition constraints 
at the end of assembly using, void ConstraintMatrix::condense.
However, it seems this function doesn't work with 
TrilinosWrappers::sparsematrix but dealii ones. 

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.

[deal.II] Periodic Boundary Condition

[Hamed Babaei]

Hi all,

Following step-45, I am going to implement periodic boundary condition on 
my code for a Thermo-elastic problem, in which the thermal and elastic 
parts are solved separately (I have two different dof_handelers for 
temprature and displacement fields without any block structure). I am 
wondering how to apply periodic constrains, made in setup_system, in the 
thermal field assemble_system for which I don't use 
"constraints.distribute_local_to_global ()" function. Thermal 
Assemble_system is the following:

  template 
  void Solid::assemble_system_eta ()
  {
system_matrix_eta = 0;
system_rhs_eta = 0;
mass_matrix = 0;
laplace_matrix = 0;
nl_matrix = 0;
nl_term = 0;

FEValues fe_values_eta (fe_eta, qf_cell_eta,
 update_values  | update_gradients |
update_quadrature_points | update_JxW_values);

FullMatrix   cell_mass_matrix(dofs_per_cell_eta, 
dofs_per_cell_eta);
FullMatrix   cell_laplace_matrix (dofs_per_cell_eta, 
dofs_per_cell_eta);
FullMatrix   cell_nl_matrix  (dofs_per_cell_eta, 
dofs_per_cell_eta);
Vector   cell_nl_term(dofs_per_cell_eta);

std::vector 
local_dof_indices(dofs_per_cell_eta);

std::vector phi (dofs_per_cell_eta);
std::vector >   grad_phi (dofs_per_cell_eta);

typename DoFHandler::active_cell_iterator
cell = dof_handler_eta.begin_active(),
endc = dof_handler_eta.end();
for (; cell!=endc; ++cell)
  if (cell->is_locally_owned())
  {
fe_values_eta.reinit(cell);

cell_mass_matrix = 0;
cell_laplace_matrix = 0;
cell_nl_matrix = 0;
cell_nl_term = 0;

PointHistory *lqph =
  
reinterpret_cast(cell->user_pointer());

   for (unsigned int q=0; qget_dof_indices (local_dof_indices);

for (unsigned int i=0; i

Re: [deal.II] Re: Convergence problem arising for large number of DoFs

[Hamed Babaei]

Hi All,

After two months of struggling with the parallel code, I finally found the 
bug. I had made a stupid mistake, initializing the temporary distributed 
solution inside the Newton loop 

I have no word to thank all of you dear friends, Wolfgang and Daniel in 
particular, for your incredible help. 

Best regards,
Hamed

-- 
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: Convergence problem arising for large number of DoFs

[Hamed Babaei]

Hi all,

It seems that although before first call to solve(in the zero newoton 
iteration) very few system matrix components are zero but after first 
solving in the first Newoton iteration, most of the system_matrix 
components are zero except all the diagonal components and few off-diagonal 
ones. Nevertheless, in the sequential code, dealii::SolverCG + SSOR can 
handle the system matrix and rhs_norm in every iteration reduces and 
Newoton solver converges. Having that said, for the parallel code, 
Trilinos::SolverCG can not converge with SSOR preconditioner but converges 
with Jacobi or AMG, however, the rhs_norm remains constant and soesn't 
shrink in every itteration so Newoton solver can not converge. The same 
story even when I used TrilinosWrappers::SolverDirect 
("Amesos_Superludist").

I was wondering why Trilinos solvers can not deal with the system_matrix 
and why Newoton iterations doesn't converge.

The system_matrix components *after* first newoton iteration is attached.

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.


System_matrix.pdf
Description: Adobe PDF document

Re: [deal.II] Re: Convergence problem arising for large number of DoFs

[Hamed Babaei]

Dear Daniel,

Then I would expect that the solver should behave the same for both 
> matrices. Are you still running into the same problems using just 4 cells 
> with your parallel code?
>

Yes, the parallel code is not solved by SolverCG+SSOR even for only 4 cell. 
It is really weird to me that despite identical system_matrices (ignoring 
numerical error) I have this problem with parallel code.!!!

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.

Re: [deal.II] Re: Convergence problem arising for large number of DoFs

[Hamed Babaei]

Dear Wolfgang,

Thanks for your help. I can reduce my problem to four 2D elements. I would 
like to plot the system_matrix of sequential and parallel codes to compare 
them but I am not sure how to do so.
Can I loop over the elements of system_matrix and simply plot them?

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.

[deal.II] Re: Convergence problem arising for large number of DoFs

[Hamed Babaei]

Hi all,

I am attaching the sequential and parallel codes as well as parameters.prm 
file. 
It would be appreciated if you could compare tangent_matrix (the 
system_matrix of the elasticity part of the problem) of parallel and 
sequential codes to find out where I have ruined the S.P.D condition that 
the parallel code is no longer solved by SolverCG+SSOR. 

The code is somehow similar to a Thermoelastic problem at large strains. 
The sequential code is based on steps 25, 18 and 44 and the parallelization 
has been done based on step-40. 

P.S. : All the matrices or vectors which have suffix "eta"  such as 
"solution_eta" or "system_matrix_eta" are related to the thermal part and 
the ones without suffix are related to elastic part.

Thank you very much for your time.

Best,
Hamed


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

#include 


namespace LA
{
#if defined(DEAL_II_WITH_PETSC) && !(defined(DEAL_II_WITH_TRILINOS) && defined(FORCE_USE_OF_TRILINOS))
  using namespace dealii::LinearAlgebraPETSc;
#  define USE_PETSC_LA
#elif defined(DEAL_II_WITH_TRILINOS)
  using namespace dealii::LinearAlgebraTrilinos;
#else
#  error DEAL_II_WITH_PETSC or DEAL_II_WITH_TRILINOS required
#endif
}

#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 PhaseField
{
  using namespace dealii;
// INPUT OF PARAMETERS
  namespace Parameters
  {
struct FESystem
{
  unsigned int poly_degree;
  unsigned int quad_order;
  static void
  declare_parameters(ParameterHandler );

  void
  parse_parameters(ParameterHandler );
};

void FESystem::declare_parameters(ParameterHandler )
{
  prm.enter_subsection("Finite element system");
  {
prm.declare_entry("Polynomial degree", "1",
  Patterns::Integer(0),
  "Displacement system polynomial order");
prm.declare_entry("Quadrature order", "2",
  Patterns::Integer(0),
  "Gauss quadrature order");
  }
  prm.leave_subsection();
}

void FESystem::parse_parameters(ParameterHandler )
{
  prm.enter_subsection("Finite element system");
  {
poly_degree = prm.get_integer("Polynomial degree");
quad_order = prm.get_integer("Quadrature order");
  }
  prm.leave_subsection();
}

struct Geometry
{
  unsigned int refinement;
  double   scale;
  static void
  declare_parameters(ParameterHandler );
  void
  parse_parameters(ParameterHandler );
};
void Geometry::declare_parameters(ParameterHandler )
{
  prm.enter_subsection("Geometry");
  {
prm.declare_entry("Global refinement", "4",
Patterns::Integer(0),
  "Global refinement level");
prm.declare_entry("Grid scale", "1e-9",
  Patterns::Double(0.0),
  "Global grid scaling factor");
  }
  prm.leave_subsection();
}
void Geometry::parse_parameters(ParameterHandler )
{
  prm.enter_subsection("Geometry");
  {
	refinement = prm.get_integer("Global refinement");
scale = prm.get_double("Grid scale");
  }
  prm.leave_subsection();
}

 /
struct Materials
{
  double lambda0; // austenite phase
  double mu0; // austenite phase
  double lambda1;  // martensite phase
  double mu1;  // martensite phase
  double L;// interface mobility
  double beta; // gradient energy coefficient
  double A0;   // parameter for barrier height
  double theta;// temperature
  double thetac;   // critical temperature
  double thetae;   // equilibrium temperature
  double thetan;   // Crank-Nicolson parameter

Re: [deal.II] Re: Convergence problem arising for large number of DoFs

[Hamed Babaei]

Hi all,

I have a quick question for you. I was wondering if the system_matrix of 
the sequential code should be exactly the same as parallel on for the same 
problem condition.
I mean if my sequential code is solved by SolverCG and SSOR preconditioner, 
the parallel one should necessarily be solved by the same solver or 
System_matrix features  may change through palatalization.

Thanks, 

On Wednesday, November 2, 2016 at 9:29:54 AM UTC-5, Bruno Turcksin wrote:
>
> Hamed, 
>
> 2016-11-02 10:22 GMT-04:00 Hamed Babaei <hamedb...@gmail.com >: 
>
> > It seems that my system_matrix is not s.p.d any more although in 
> sequential 
> > code it was. 
> Then, there is a bug in your code. You can compare the matrix produced 
> by the sequential code with the one in parallel do help you find the 
> bug. 
>
> 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.

Re: [deal.II] Re: Convergence problem arising for large number of DoFs

[Hamed Babaei]

Dear Wolfgang,

I wrote the sequential code first and It converged to the exact solution 
using SolverCG and SSOR preconditioner. 
However, the parallel code diverges at the first call to SolverCG but using 
GMRES or Bicgstab it works though up to 15 DoFs not more.

First and the most important question to me is what parallelization does to 
my system_matrix that it is not solved by SolverCG anymore. I guess the 
source of the problem is there.
Second, why it is solved by GMRES and Bicgstab for few DoFs (up to 15) 
and not more( for example, 1 million). 

It seems that my system_matrix is not s.p.d any more although in sequential 
code it was. 

Thanks for your time,
Hamed

On Wednesday, November 2, 2016 at 8:00:23 AM UTC-5, Wolfgang Bangerth wrote:
>
> On 11/01/2016 10:31 AM, Hamed Babaei wrote: 
> > 
> > I mean non of the standard preconditioners helped to make it converged. 
>
> Like Bruno mentioned, if the standard preconditioners don't work, then you 
> will need to think about more complicated ones. Take a look at my video 
> lectures on preconditioners, for example. 
>
> There is of course also the possibility that your matrix is not correct 
> and 
> indeed singular. In that case, you *have* to expect problems. Have you 
> verified that your numerical solutions converge to the exact ones for 
> small 
> problems? 
>
> 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.

Re: [deal.II] Re: Convergence problem arising for large number of DoFs

[Hamed Babaei]

Dear Bruno,

I mean non of the standard preconditioners helped to make it converged.

On Monday, October 31, 2016 at 3:13:35 PM UTC-5, Bruno Turcksin wrote:
>
> Hamed, 
>
> 2016-10-31 15:56 GMT-04:00 Hamed Babaei <hamedb...@gmail.com >: 
>
> > I have already checked SSOR and AMG non of them make a difference . ILU 
> > doesn't work and gives this error: "[0]PETSC ERROR: MatSolverPackage 
> petsc 
> > does not support matrix type mpiaij" 
> What do you mean they don't make a difference? You should get a 
> different residual after 1 steps for each preconditioner. 
>
> > I don't know how to increase number of iteration. Although I have used 
> the 
> > "SolverControl solver_control (dof_handler.n_dofs(), 
> > 1e-12*system_rhs.l2_norm());"  as solver control and my n_dof is much 
> more 
> > than 1 but it doesn't go more than 1 itteration !!! 
> I don't know. This could be a limit in PETSc? I never use PETSc so I don't 
> know. 
>
> > I was wondering which other preconditioner to use and how to increase 
> number 
> > of iterations . 
> If standard preconditioners don't work very well for your problem, you 
> need a physics-based preconditioner. You will need to find and 
> implement it yourself since it needs to be tailored for your equation. 
> Something that you can also try, it's to restart less often GMRES 
> (
> http://www.dealii.org/developer/doxygen/deal.II/structPETScWrappers_1_1SolverGMRES_1_1AdditionalData.html)
>  
>
> Restarting less often will make GMRES converge faster but you will 
> need more memory. 
>
> 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: Convergence problem arising for large number of DoFs

[Hamed Babaei]

Dear Bruno,

I have already checked SSOR and AMG non of them make a difference . ILU 
doesn't work and gives this error: "[0]PETSC ERROR: MatSolverPackage petsc 
does not support matrix type mpiaij"

I don't know how to increase number of iteration. Although I have used the 
"SolverControl solver_control (dof_handler.n_dofs(), 
1e-12*system_rhs.l2_norm());"  as solver control and my n_dof is much more 
than 1 but it doesn't go more than 1 itteration !!!

I was wondering which other preconditioner to use and how to increase 
number of iterations .

Thanks



On Monday, October 31, 2016 at 2:14:38 PM UTC-5, Bruno Turcksin wrote:
>
> Hamed,
>
> On Monday, October 31, 2016 at 3:01:20 PM UTC-4, Hamed Babaei wrote:
>>
>> This error message can indicate that you have simply not allowed a 
>> sufficiently large number of iterations for your iterative solver to 
>> converge. This often happens when you increase the size of your problem. In 
>> such cases, the last residual will likely still be very small, and you can 
>> make the error go away by increasing the allowed number of iterations when 
>> setting up the SolverControl object that determines the maximal number of 
>> iterations you allow.
>>
>
> BlockJacobi is a bad preconditioner, you will need increase the maximal 
> number of iterations and/or use a better preconditioner.
>
> 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] Convergence problem arising for large number of DoFs

[Hamed Babaei]

Hi all,

I have parallelized a code in which I use GMRES or Bicgstab solvers with 
BlockJacobi preconditioner for a symmetric but not possitive-difinite 
system matrix. The problem is that it converges for 15 DoFs and fewer 
(on 1 node with 16 processor) but does'nt converge for larger number of 
DoFs for example 50 or one million (even with 64 processors). 

It would be appreciated if you could give me any clue to find out how to 
make it converge for large number of DoFs. The following is the error I 
usually get.

Thanks 
Hamed


Exception on processing:


An error occurred in line <151> of file 

 
in function
void dealii::PETScWrappers::SolverBase::solve(const 
dealii::PETScWrappers::MatrixBase&, dealii::PETScWrappers::VectorBase&, 
const dealii::PETScWrappers::VectorBase&, const 
dealii::PETScWrappers::PreconditionerBase&)
The violated condition was:
false
The name and call sequence of the exception was:
SolverControl::NoConvergence (solver_control.last_step(), 
solver_control.last_value())
Additional Information:
Iterative method reported convergence failure in step 1. The residual 
in the last step was 1.54204e-06.

This error message can indicate that you have simply not allowed a 
sufficiently large number of iterations for your iterative solver to 
converge. This often happens when you increase the size of your problem. In 
such cases, the last residual will likely still be very small, and you can 
make the error go away by increasing the allowed number of iterations when 
setting up the SolverControl object that determines the maximal number of 
iterations you allow.

The other situation where this error may occur is when your matrix is not 
invertible (e.g., your matrix has a null-space), or if you try to apply the 
wrong solver to a matrix (e.g., using CG for a matrix that is not symmetric 
or not positive definite). In these cases, the residual in the last 
iteration is likely going to be large.


Aborting!

ERROR: Uncaught exception in MPI_InitFinalize on proc 1. Skipping 
MPI_Finalize() to avoid a deadlock.

-- 
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: what do SolverControl do for Direct Solvers?

[Hamed Babaei]

Dear Bruno,

If I solve a problem with one million DOFs on a cluster of 64 processors, 
every processor will be dealing with around 15000 DOFs. 
So you mean even in that case direct solver would be slow? I thought direct 
solver would solve the problem only on those DOFs that it owns.

Best,
Hamed

On Friday, October 28, 2016 at 7:59:06 AM UTC-5, Bruno Turcksin wrote:
>
> Hamed, 
>
> 2016-10-27 19:03 GMT-04:00 Hamed Babaei <hamedb...@gmail.com >: 
>
> > I've been using the default one, Amesos_KLU, and I checked 
> > Amesos_Superludist but trilinos doesn't recognize it. 
> > when we installed dealii and its dependencies on the cluster, I think we 
> > didn't uncomment the "#PACKAGES="${PACKAGES} once:superlu_dist" in 
> candi.cfg 
> > so it seems we need to install them again. 
> Yes, Superlu_dist is a code independent of Trilinos. You will need to 
> install it and reinstall Trilinos 
>
> > I want to solve my problem with more than one million DOFs, So you mean 
> > Amesos_Superludist will work for that many DOFs when running on cluster 
> ? 
> It will work but it will take a while. With that many DOFs, you should 
> use an iterative solver. Superlu_dist is a lot faster than KLU and it 
> will be good for the current size of your problem but when you have 
> more than one million of DOFs it will be slow. 
>
> 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: what do SolverControl do for Direct Solvers?

[Hamed Babaei]

Dear Bruno,

I've been using the default one, Amesos_KLU, and I checked 
Amesos_Superludist but trilinos doesn't recognize it.
when we installed dealii and its dependencies on the cluster, I think we 
didn't uncomment the "#PACKAGES="${PACKAGES} once:superlu_dist" in 
candi.cfg so it seems we need to install them again.

I want to solve my problem with more than one million DOFs, So you mean 
Amesos_Superludist will work for that many DOFs when running on cluster ?

Thanks


On Thursday, October 27, 2016 at 5:37:31 PM UTC-5, Bruno Turcksin wrote:
>
> Hamed,
>
> On Thursday, October 27, 2016 at 5:40:04 PM UTC-4, Hamed Babaei wrote:
>
>> First: I have parallelized a code in which I use 
>> TrilinosWrappers::SolverDirect 
>> . The problem is that with around 2 DOFs the code is run on my own 
>> machine and on the cluster as well but when I increase DOF even to 15 
>> it get stuck at the first call to direct solver without giving any error 
>> and Direct solver never got done. 
>>
> 15 is a lot of unknowns for a direct solver. 10 unknowns is around 
> the time when you need to switch to an iterative method to get the results 
> quickly enough. Which solver are you using? If you are using Amesos_KLU, 
> the solver is serial so using more processors is pointless. You need to use 
> Amesos_Superludist.
>  
>
>> Second:  I don't know what solver control do for SolverDirect !!  
>> SolverDirect 
>> <https://www.dealii.org/8.4.1/doxygen/deal.II/classTrilinosWrappers_1_1SolverDirect.html#a5d216dac60a7bdf696c206d926974d1f>
>>  (SolverControl 
>> <https://www.dealii.org/8.4.1/doxygen/deal.II/classSolverControl.html> , 
>> const AdditionalData 
>> <https://www.dealii.org/8.4.1/doxygen/deal.II/structTrilinosWrappers_1_1SolverDirect_1_1AdditionalData.html>
>>  =AdditionalData 
>> <https://www.dealii.org/8.4.1/doxygen/deal.II/structTrilinosWrappers_1_1SolverDirect_1_1AdditionalData.html>
>> ()) 
>> I thought SolverControl has meaning when we have iterative solvers.
>>
> You are right. It is use internally to check that the direct solver was 
> able to find the solution but the parameters that you put in SolverControl 
> are not used. 
>
>>
>> I was wondering why the code doesn't work with increased DOF and if I can 
>> do something to find out where the problem is and how to solve it.
>>
> It probably works but you may need to wait a long time for the direct 
> solver to solve the system.
>
> 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: what do SolverControl do for Direct Solvers?

[Hamed Babaei]

I want to add that for running on Cluster I used up to 64 processors (4 
nodes, each with 16 processor).  

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

[deal.II] what do SolverControl do for Direct Solvers?

[Hamed Babaei]

Hi all,

I have two question which may be somehow related:

First: I have parallelized a code in which I use TrilinosWrappers::SolverDirect 
. The problem is that with around 2 DOFs the code is run on my own 
machine and on the cluster as well but when I increase DOF even to 15 
it get stuck at the first call to direct solver without giving any error 
and Direct solver never got done. 
Second:  I don't know what solver control do for SolverDirect !!  
SolverDirect 

 (SolverControl 
 , 
const AdditionalData 

 =AdditionalData 

()) 
I thought SolverControl has meaning when we have iterative solvers.

I was wondering why the code doesn't work with increased DOF and if I can 
do something to find out where the problem is and how to solve it.

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.

Re: [deal.II] How to addapt MappingQEulerian class for Parallel codes using the PETSc wrapper classes

[Hamed Babaei]

Dear Timo and Wolfgang,

Thank you very much for your guidance.

I should have passed the ghosted solution vector itself as the third 
parameter in MappingQEulerian instead of a non-ghosted copy of it.

Best, Regards,
Hamed

On Monday, October 24, 2016 at 3:37:48 PM UTC-5, Timo Heister wrote:
>
> > It would be appreciated if you could give me any clue to resolve the 
> error. 
>
> You need a ghosted vector as the parameter in MappingQEulerian. Is 
> that the case for you? 
>
> See 
> https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_dealii_dealii_blob_master_tests_mappings_mapping-5Fq-5Feulerian-5F07.cc=CwIBaQ=Ngd-ta5yRYsqeUsEDgxhcqsYYY1Xs5ogLxWPA_2Wlc4=4k7iKXbjGC8LfYxVJJXiaYVu6FRWmEjX38S7JmlS9Vw=yBGGt1DGkK9Yl7pJR_GSSnJwOsJafwPvtA_WoS_0lj0=vZ5fLIcJ2VyttPXyTy7O68tEVbMToLqcwyhrli5jmtE=
>  
> for a simple example. 
>
> -- 
> Timo Heister 
> http://www.math.clemson.edu/~heister/ 
>

-- 
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] Direct Solvers in Parallel

[Hamed Babaei]

Hi friends,

I am parallelizing a code similar to step-44 in which it is possible to use 
either an iterative solver SolverCG or a direct solver, 
SparseDirectUMFPACK. I have used the latter in the non-parallel code and 
works great.
Using iterative solvers like SolverCG I have problem in convergence so I 
want to check a direct solver which works in parallel. My problem is that 
my code doesn't recognize  PETScWrappers::SparseDirectMUMPS 

  
nor TrilinosWrappers::SolverDirect 

 .
I have installed Dealii and all of its dependent libraries (Petsc, 
Trilinos, P4est ...) via Candi (https://github.com/koecher/candi). I was 
wondering which direct solver I should use which works the same as 
SparseDirectUMFPACK and how to make dealii know them.

Thanks,
Hamed

-- 
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 addapt MappingQEulerian class for Parallel codes using the PETSc wrapper classes

[Hamed Babaei]

Dear Timo,

Regarding using MappingQEulerian class in parallel, I did what you said for 
the vector type template argument. It works when I run the code on my own 
machine but I got the following error when I run it on cluster:

An error occurred in line <1480> of file 
 
in function
dealii::IndexSet::size_type 
dealii::IndexSet::index_within_set(dealii::IndexSet::size_type) const
The violated condition was:
is_element(n) == true
The name and call sequence of the exception was:
ExcIndexNotPresent (n)
Additional Information:
The global index 0 is not an element of this set.

Stacktrace:
---
#0  ./parallel_large-strain_PT: dealii::IndexSet::index_within_set(unsigned 
int) const
#1  /shared/hpc/deal.ii-candi/deal.II-v8.4.1/lib/libdeal_II.g.so.8.4.1: 
dealii::PETScWrappers::internal::VectorReference::operator double() const
#2  ./parallel_large-strain_PT: PhaseField::Solid<3>::output_results() const
#3  ./parallel_large-strain_PT: PhaseField::Solid<3>::run()
#4  ./parallel_large-strain_PT: main

Following is my out_put object. problem is related to the MappingQEulerian 
since when I uncomment it (the red lines) I don't get error.

 template 
  void Solid::output_results() const
  {
DataOut data_out;

std::vector displacement_names;
switch (dim)
  {
  case 1:
displacement_names.push_back ("displacement");
  break;
  case 2:
  displacement_names.push_back ("x_displacement");
  displacement_names.push_back ("y_displacement");
  break;
  case 3:
  displacement_names.push_back ("x_displacement");
  displacement_names.push_back ("y_displacement");
  displacement_names.push_back ("z_displacement");
  break;
 default:
  Assert (false, ExcNotImplemented());
 }

data_out.add_data_vector (dof_handler, solution, displacement_names);
 data_out.add_data_vector (dof_handler_eta, solution_eta, 
"order_parameter");
LA::MPI::Vector soln ;
soln.reinit(locally_owned_dofs,mpi_communicator);
for (unsigned int i = 0; i < soln.size(); ++i)
  soln(i) = solution(i);
MappingQEulerian q_mapping(degree, dof_handler, 
soln);

Vector subdomain (triangulation.n_active_cells());
for (unsigned int i=0; i
> > "To enable the use of the MappingQ1Eulerian class also in the context of 
> > parallel codes using the PETSc wrapper classes, the type of the vector 
> can 
> > be specified as template parameter EulerVectorType Not specifying this 
> > template argument in applications using the PETSc vector classes leads 
> to 
> > the construction of a copy of the vector which is not accessible 
> > afterwards!" 
> > I can't understatnd what it means and what changes I should make. 
>
> The class MappingQ1Eulerian has a template argument for the vector 
> type. You need to set this to be a parallel vector (like 
> PETScWrappers::MPI::Vector). 
>
> -- 
> Timo Heister 
> http://www.math.clemson.edu/~heister/ 
>

-- 
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 apply boundary values for a particular point on the boundary instead of the whole boundary surface

[Hamed Babaei]

Dear Stephen,
Thank you very much for your incredible help.

On Saturday, October 22, 2016 at 6:43:29 AM UTC-5, Stephen DeWitt wrote:
>
> In case that link dies, here's the snippet of code:
>
>
> // Set constraints to pin the solution if there are no Dirichlet BCs for a 
> component of a variable in an elliptic equation
> template 
> void MatrixFreePDE::setRigidBodyModeConstraints( std::vector 
> rigidBodyModeComponents, ConstraintMatrix * constraints, DoFHandler* 
> dof_handler){
>
> if ( rigidBodyModeComponents.size() > 0 ){
>
>   // Choose the point where the constraint will be placed. Must be the 
> coordinates of a vertex.
>   dealii::Point target_point(0,0);
>
>   unsigned int vertices_per_cell=GeometryInfo::vertices_per_cell;
>
>   // Loop over each locally owned cell
>   typename DoFHandler::active_cell_iterator cell= 
> dof_handler->begin_active(), endc = dof_handler->end();
>
>   for (; cell!=endc; ++cell){
>   if (cell->is_locally_owned()){
>   for (unsigned int i=0; i<vertices_per_cell; ++i){
>
>   // Check if the vertex is the target vertex
>   if (target_point.distance (cell->vertex(i)) < 
> 1e-2 * cell->diameter()){
>
>   // Loop through the list of components with 
> rigid body modes and add an inhomogeneous constraint for each
>   for (unsigned int component_num = 0; 
> component_num < rigidBodyModeComponents.size(); component_num++){
>   unsigned int 
> nodeID=cell->vertex_dof_index(i,component_num);
>   constraints->add_line(nodeID);
>   
> constraints->set_inhomogeneity(nodeID,0.0);
>   }
>  }
>  }
>  }
>}
> }
> }
>
>
>
>
>
> On Saturday, October 22, 2016 at 6:49:41 AM UTC-4, Stephen DeWitt wrote:
>>
>> Hi Hamed,
>> I was in the same boat a few weeks ago. After reading through some of the 
>> other posts on this list (which Jean-Paul linked), I wrote the following 
>> method that you and maybe others will find useful.
>>
>> It takes in a vector containing the list of components of a field that 
>> need point constraints (rigidBodyModeComponents) and adds a constraint 
>> where needed:
>>
>> https://github.com/prisms-center/phaseField/blob/next/src/matrixfree/boundaryConditions.cc
>> (lines 43-75)
>>
>> Cheers!
>> Steve
>>
>>
>> On Thursday, October 20, 2016 at 7:38:19 PM UTC-4, Hamed Babaei wrote:
>>>
>>> Hi friends,
>>>
>>> For an elastic problem, I am going to apply zero boundary displacements 
>>> for three specific points on the center of -x, -y and -z planes of a cubic 
>>> domain.
>>> I have already done this but for the boundary surface not a boundary 
>>> point (the same as incremental_boundary_displacement in step-18). 
>>> The following is what I wrote to do so which doesn't work properly. In 
>>> fact, it fixes the whole -x, -y and -z surfaces, not just for the three 
>>> points on them that I intended. 
>>>
>>>   template 
>>>   class BoundaryCondition :  public Function
>>>   {
>>>   public:
>>>  BoundaryCondition (const int boundary_id);
>>> virtual void vector_value (const Point ,
>>> Vector   ) const;
>>> virtual void vector_value_list (const std::vector<Point 
>>> > ,
>>>std::vector<Vector >   
>>> _list) const;
>>>   private:
>>> const int boundary_id;
>>>
>>>   };
>>>   template 
>>>   BoundaryCondition::BoundaryCondition (const int 
>>> boundery_id)
>>> :
>>> Function (dim),
>>> boundary_id(boundary_id)
>>>
>>>   {}
>>>   template 
>>>   inline
>>>   void
>>> BoundaryCondition::vector_value (const Point ,
>>> Vector   ) const
>>>   {
>>> Assert (values.size() == dim,
>>> ExcDimensionMismatch (values.size(), dim));
>>>
>>> Point point_x;
>&

Re: [deal.II] Re: How to use Trilinos instead of Petsc in step-40

[Hamed Babaei]

Dear Daniel,

Please look at the pdf file attached. Thanks

On Friday, October 21, 2016 at 4:51:22 PM UTC-5, Daniel Arndt wrote:
>
> Hamed,
>
> At least I can't open your file in a way to clearly see what your 
> bilinearform is, but could just guess from the implementation. Can you 
> either write it here (using LaTeX or similar) or upload the file as pdf 
> again?
>
> 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.


Bilinearform.pdf
Description: Adobe PDF document

Re: [deal.II] Re: How to use Trilinos instead of Petsc in step-40

[Hamed Babaei]

Dear Wolfgang,

In attach you shall find the bilinear form and the assemble part of the 
code as well.

Thanks,

On Friday, October 21, 2016 at 11:14:52 AM UTC-5, Wolfgang Bangerth wrote:
>
>
> > On Thursday, October 20, 2016 at 9:48:55 PM UTC-5, Wolfgang Bangerth 
> wrote: 
> > 
> > What is the norm of the first right hand side? You now set the 
> tolerance in 
> > 
> > the first iteration to a fixed value of 1e-16, but how does this 
> compare to 
> > the previous value of 1e-12*system_rhs.l2_norm()? 
> > 
> > 
> > The system_rhs.l2_norm() before the first call to CGSolverat in the zero 
> time 
> > step was 1.82492e-09, so that the 1e-12*system_rhs.l2_norm() became in 
> the 
> > order of 1e-21 by which CGSOlver didn't converge. 
>
> OK. I have no suggestion, though, why one or the other would lead to 
> convergence or non-convergence of CG. 
>
>
> > That still seems wrong to me. As I mentioned in a previous email, 
> you 
> > ought to 
> > make sure that the matrix you build is the matrix you *want* to 
> build. 
> > 
> > 
> > I was wondering how to make sure the matrix I built is what I meant. 
>
> You need to compare properties of your bilinear form with those of your 
> matrix. For example, if your bilinear form is symmetric, then your matrix 
> is 
> as well -- and that is something you can test. If your bilinear form is 
> positive definite, then your matrix needs to be as well and as a 
> consequence 
> CG needs to converge. If CG does not seem to converge but your bilinear 
> form 
> is s.p.d., you know you have a bug in your matrix assembly. If CG 
> converges on 
> one processor but not on two for the same problem, then you know you have 
> a 
> bug in the parallel matrix assembly. Etc. 
>
> Since you haven't stated your bilinear form or the problem you are 
> solving, I 
> can't tell you whether the matrix is supposed to be s.p.d. 
>
> 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.


Bilinearform.docx
Description: MS-Word 2007 document

Re: [deal.II] Re: How to use Trilinos instead of Petsc in step-40

[Hamed Babaei]

Dear Wolfgang,

On Thursday, October 20, 2016 at 9:48:55 PM UTC-5, Wolfgang Bangerth wrote:
>
> What is the norm of the first right hand side? You now set the tolerance 
> in 
>
the first iteration to a fixed value of 1e-16, but how does this compare to 
> the previous value of 1e-12*system_rhs.l2_norm()? 
>

The system_rhs.l2_norm() before the first call to CGSolverat in the zero 
time step was 1.82492e-09, so that the 1e-12*system_rhs.l2_norm() became in 
the order of 1e-21 by which CGSOlver didn't converge.

>
> That still seems wrong to me. As I mentioned in a previous email, you 
> ought to 
> make sure that the matrix you build is the matrix you *want* to build. 
>
 
I was wondering how to make sure the matrix I built is what I meant. 
In addition, Bruno already mentioned sometimes PreconditionerAMG has 
problem with non-positive definite System_matrices. 

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.

[deal.II] How to apply boundary values for a particular point on the boundary instead of the whole boundary surface

[Hamed Babaei]

Hi friends,

For an elastic problem, I am going to apply zero boundary displacements for 
three specific points on the center of -x, -y and -z planes of a cubic 
domain.
I have already done this but for the boundary surface not a boundary point 
(the same as incremental_boundary_displacement in step-18). 
The following is what I wrote to do so which doesn't work properly. In 
fact, it fixes the whole -x, -y and -z surfaces, not just for the three 
points on them that I intended. 

  template 
  class BoundaryCondition :  public Function
  {
  public:
 BoundaryCondition (const int boundary_id);
virtual void vector_value (const Point ,
Vector   ) const;
virtual void vector_value_list (const std::vector 
,
   std::vector   
_list) const;
  private:
const int boundary_id;

  };
  template 
  BoundaryCondition::BoundaryCondition (const int boundery_id)
:
Function (dim),
boundary_id(boundary_id)

  {}
  template 
  inline
  void
BoundaryCondition::vector_value (const Point ,
Vector   ) const
  {
Assert (values.size() == dim,
ExcDimensionMismatch (values.size(), dim));

Point point_x;
point_x(1) = 5;
point_x(2) = 5;

Point point_y;
point_y(0) = 5;
point_y(2) = 5;

Point point_z;
point_z(0) = 5;
point_z(1) = 5;


if  (boundary_id ==0 && ((p-point_x).norm_square() 
<(0.5e-9)*(0.5e-9)))
   values(0) = 0;
else if (boundary_id ==2 && ((p-point_y).norm_square() < 
(0.5e-9)*(0.5e-9)))
   values(1)= 0;
else if (boundary_id ==4 && ((p-point_z).norm_square() < 
(0.5e-9)*(0.5e-9)))
   values(2)= 0;

  }
  template 
  void
BoundaryCondition::vector_value_list (const 
std::vector ,
 std::vector   
_list) const
  {
const unsigned int n_points = points.size();
Assert (value_list.size() == n_points,
ExcDimensionMismatch (value_list.size(), n_points));
for (unsigned int p=0; p

Re: [deal.II] Re: How to use Trilinos instead of Petsc in step-40

[Hamed Babaei]

Hi all,

I just made some changes in the second parameter of SolverControl in the 
following way which apparently resolved the problem:

my problem is nonlinear so I solve it by newton method. It seems that in 
the first Newton-Rophson iteration the CG solver can not converge if the 
tolerance for success of CG iterations is considered 
"1e-12*system_rhs.l2_norm()" 
but It converges with a smaller number like 10e-16. what I did to 
solver_control follows:

double tol ;
 if (newton_rophson_itteration_number==0)
   tol = 1e-16;
else
tol =1e-12*system_rhs.l2_norm();
 
SolverControl solver_control (dof_handler.n_dofs(), tol);

and the following is system_rhs.l2_norm (newton residual) for the first 
time step:

newton iteration   system_rhs.l2_normnumber of CG iterations
0 1.82492e-09 97
1 7.71201e-07 185
2 2.1306e-07   207
3 6.89427e-08 199
4 3.16879e-09 198
5 8.83219e-12 201
6 9.55033e-17 202

As you can see in the first newton iteration residual increases then starts 
to decrease till convergence criteria is satisfied. 
Although the problem has been resolved by this trick but I don't know why I 
need to consider different convergence criteria for CG solver in the first 
newton iteration.

In addition, the CG solver still doesn't converge with PreconditionAMG but 
it works with PreconditionJacobi. 

Thanks,

On Wednesday, October 19, 2016 at 11:01:05 AM UTC-5, Bruno Turcksin wrote:
>
> 2016-10-19 11:51 GMT-04:00 Hamed Babaei <hamedb...@gmail.com >: 
>
> > I need to compute the determinant of system_matrix to check if its 
> positive 
> > or negative so that I know my system_matrix is positive definite or not. 
> That's what I thought. You don't need to know the determinant to do 
> that, you can check the sign of the eigenvalues. If all the 
> eigenvalues are positive, the matrix is positive definite. You can 
> either use a separate routine to compute the eigenvalues (see step-36 
> dealii.org/8.4.1/doxygen/deal.II/step_36.html) or you can get an 
> estimate of the eigenvalues when you are using GMRES (see 
>
> http://dealii.org/8.4.1/doxygen/deal.II/classSolverGMRES.html#a7b12c13beeb2d307f244e4dd48c06859)
>  
>
>
> 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.

Re: [deal.II] Re: How to use Trilinos instead of Petsc in step-40

[Hamed Babaei]

I need to compute the determinant of system_matrix to check if its positive 
or negative so that I know my system_matrix is positive definite or not.

On Wednesday, October 19, 2016 at 10:38:09 AM UTC-5, Bruno Turcksin wrote:
>
> Hamed,
>
> On Wednesday, October 19, 2016 at 10:25:37 AM UTC-4, Hamed Babaei wrote:It 
> seems that I've been choosing wrong solver for my problem. Since my problem 
> has some sort of instability in its nature, 
>>
>> although my system_matrix is symmetric, it is not necessarily positive 
>> definite. So I have two questions: 
>> First, Is there any function in Dealii or Petsc to compute the 
>> determinant of system_matrix. I have used a determinant function for 3*3 
>> matrices and it doesn't work for higher dimensions.
>>
> Why do you need this?
>  
>
>> Second, which solver and preconditioner can be used for symmetric but 
>> non-positive definite system_matrix.
>>
> For the solver, you can use GMRES, BiCGSTAB, or MinRES. For the 
> preconditioner, you can use ILU, SSOR, or algebraic multigrid (AMG). 
> Sometimes AMG have difficulties with non-s.p.d. matrices though.
>
> 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.

Re: [deal.II] Re: How to use Trilinos instead of Petsc in step-40

[Hamed Babaei]

Hi everyone,

Thank you all for your incredible guides.
It seems that I've been choosing wrong solver for my problem. Since my 
problem has some sort of instability in its nature, 
although my system_matrix is symmetric, it is not necessarily positive 
definite. So I have two questions: 
First, Is there any function in Dealii or Petsc to compute the determinant 
of system_matrix. I have used a determinant function for 3*3 matrices and 
it doesn't work for higher dimensions.
Second, which solver and preconditioner can be used for symmetric but 
non-positive definite system_matrix.

Thanks


On Tuesday, October 18, 2016 at 9:49:03 AM UTC-5, Wolfgang Bangerth wrote:
>
> > which solvers do you think I can examin and see if they can be used 
> > instead of SolverCG? 
>
> That's the wrong question. If you know that the matrix you have *should 
> be* 
> symmetric and positive definite, but the CG solver does not converge, then 
> you 
> need to debug why the matrix *is not* s.p.d. Using a different solver 
> might 
> allow you to solve the linear system, but since the linear system is 
> apparently not the correct one, the solution will also not be correct. 
>
> 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.

[deal.II] Re: How to use Trilinos instead of Petsc in step-40

[Hamed Babaei]

Another point I need to mention is that I had the same problem when as my 
first try of parallelization, I paralleled a much simpler code, the 
step-25, based on the step-40 . 
I received the same Segmentation Violation error from Petsc. At that time, 
I replaced PreconditionerAMG with PreconditionerJacobi and the problem 
resolved. 

However, Now even using PreconditionerJacobi can not remedy the issue even 
if it works for very few initial value and boundary conditions.

On Sunday, October 16, 2016 at 10:39:49 AM UTC-5, Daniel Arndt wrote:
>
> In case there is an error with Trilinos, what is it?
>
> Am Sonntag, 16. Oktober 2016 17:27:26 UTC+2 schrieb Daniel Arndt:
>>
>> Hamed,
>>
>>
> It still works with Petsc and Trilinos is not implemented at all. I 
> appreciate it if you could let me know what changes should be made in 
> step-40 to make it work with Trilinos.
>
 What exactly, do you mean by that? So you have both PETSc and Trilinos 
>> installed and step-40 works. Are you saying that PETSc objects are used 
>> even if you use 
>> "using namespace ::LinearAlgebraTrilinos;" ? How did you find out?
>>
>> 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] Problem with PETScWrappers::SolverCG

[Hamed Babaei]

Hi freinds,

I have parallelized a code for a Thermoelastic problem based on step-40 . 
Although the non-parallel code works for all the initial values and 
boundary conditions, the parallel one fails giving this error :
"PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, 
probably memory access out of range"

It seems from debugger backtracethat that error comes from the 
PETScWrappers::SolverCG. 
Surprisingly, the problem converges for few initial value or boundary 
conditions (not for all problem conditions) when I change the 
preconditioner from PreconditionAMG to PreconditionJacobi.

It would be appreciated if you could give me any clue why I have this 
error. The following is my solve object:

template 
 unsigned int
 Solid::solve ()
{

  LA::MPI::Vector
  completely_distributed_solution (locally_owned_dofs, 
mpi_communicator);

  SolverControl solver_control (dof_handler.n_dofs(), 1e-12);

  #ifdef USE_PETSC_LA
  LA::SolverCG solver(solver_control, mpi_communicator);
  #else
  LA::SolverCG solver(solver_control);
  #endif

  LA::MPI::PreconditionAMG preconditioner;
  LA::MPI::PreconditionAMG::AdditionalData data;
  #ifdef USE_PETSC_LA
  data.symmetric_operator = true;
  //data.strong_threshold = 0.5;
  #else
  /* Trilinos defaults are good */
  #endif
   
preconditioner.initialize(tangent_matrix, data);
   

//  LA::MPI::PreconditionJacobi preconditioner(tangent_matrix);


  solver.solve (tangent_matrix, completely_distributed_solution, 
system_rhs,
preconditioner);

  constraints.distribute (completely_distributed_solution);

  solution_update = completely_distributed_solution;

  return solver_control.last_step();
}

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.

[deal.II] Re: How to use Trilinos instead of Petsc in step-40

[Hamed Babaei]

Hi Jean-Paul,

I have installed latest version of Dealii (8.4.1) and all of its library 
dependencies (Petsc, P4est, Trilinos,...) all together via Candi ( the bash 
script installer tool) https://github.com/koecher/candi. 

Thanks

On Saturday, October 15, 2016 at 12:25:44 AM UTC-5, Jean-Paul Pelteret 
wrote:
>
> Hi Hamad,
>
> Can you please tell us which version of deal.II you are using, and provide 
> us with the installation logs?
>
> Regards,
> J-P
>
> On Friday, October 14, 2016 at 11:51:09 PM UTC+2, Hamed Babaei wrote:
>>
>> Hi friends,
>>
>> I want to use Trilinos instead of Petsc in my parallel code which is 
>> based on step-40. I guess some changes should be made in the following 
>> lines of the include files part:
>>
>> namespace LA
>> {
>> #if defined(DEAL_II_WITH_PETSC) && !(defined(DEAL_II_WITH_TRILINOS) && 
>> defined(FORCE_USE_OF_TRILINOS))
>> using namespace ::LinearAlgebraPETSc 
>> <http://dealii.org/8.4.1/doxygen/deal.II/namespaceLinearAlgebraPETSc.html>
>> ;
>> # define USE_PETSC_LA
>> #elif defined(DEAL_II_WITH_TRILINOS)
>> using namespace ::LinearAlgebraTrilinos 
>> <http://dealii.org/8.4.1/doxygen/deal.II/namespaceLinearAlgebraTrilinos.html>
>> ;
>> #else
>> # error DEAL_II_WITH_PETSC or DEAL_II_WITH_TRILINOS required
>> #endif
>> }
>>
>> My problem is even if I replace the above lines with the following:
>>
>> namespace LA
>> {
>> using namespace ::LinearAlgebraTrilinos 
>> <http://dealii.org/8.4.1/doxygen/deal.II/namespaceLinearAlgebraTrilinos.html>
>> ;
>> }
>>
>> It still works with Petsc and Trilinos is not implemented at all. I 
>> appreciate it if you could let me know what changes should be made in 
>> step-40 to make it work with Trilinos.
>>
>> Thanks,
>> Hamed
>>
>

-- 
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 use Trilinos instead of Petsc in step-40

[Hamed Babaei]

Hi friends,

I want to use Trilinos instead of Petsc in my parallel code which is based 
on step-40. I guess some changes should be made in the following lines of 
the include files part:

namespace LA
{
#if defined(DEAL_II_WITH_PETSC) && !(defined(DEAL_II_WITH_TRILINOS) && 
defined(FORCE_USE_OF_TRILINOS))
using namespace ::LinearAlgebraPETSc 
;
# define USE_PETSC_LA
#elif defined(DEAL_II_WITH_TRILINOS)
using namespace ::LinearAlgebraTrilinos 

;
#else
# error DEAL_II_WITH_PETSC or DEAL_II_WITH_TRILINOS required
#endif
}

My problem is even if I replace the above lines with the following:

namespace LA
{
using namespace ::LinearAlgebraTrilinos 

;
}

It still works with Petsc and Trilinos is not implemented at all. I 
appreciate it if you could let me know what changes should be made in 
step-40 to make it work with Trilinos.

Thanks,
Hamed

-- 
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] Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range

[Hamed Babaei]

Hi Daniel,

I hope you are doing well. I've been struggling to run step-40 using 
Trilinos instead of Petsc. 
I appreciate it if you could let me know which parts of the code should be 
changed and how.

Documentation in the include part says : " uncomment the following #define 
if you have PETSc and Trilinos installed and you prefer using Trilinos in 
this example: #define FORCE_USE_OF_TRILINOS"
I can't understand what it means by the change in the following part:

namespace LA
{
#if defined(DEAL_II_WITH_PETSC) && !(defined(DEAL_II_WITH_TRILINOS) && 
defined(FORCE_USE_OF_TRILINOS))
using namespace ::LinearAlgebraPETSc 
;
# define USE_PETSC_LA
#elif defined(DEAL_II_WITH_TRILINOS)
using namespace ::LinearAlgebraTrilinos 

;
#else
# error DEAL_II_WITH_PETSC or DEAL_II_WITH_TRILINOS required
#endif
}

Thanks

On Monday, July 4, 2016 at 5:00:51 PM UTC-5, Daniel Arndt wrote:
>
> Ehsan,
>
> All I can say: After switching the order of arguments in 
> SparseMatrix::add, your code runs for me with a recent developer version 
> and Trilinos at least.
>
> Best,
> Daniel
>
> Am Montag, 4. Juli 2016 18:59:05 UTC+2 schrieb Ehsan Esfahani:
>>
>> Dear Professor Bangerth,
>>
>> Thanks for your response. Yes, I did. As I mentioned, I got a backtrace 
>> in the debugger (eclipse) and I find out that the problem is in the line I 
>> have mentioned but I couldn't find out what's the problem in that line of 
>> the code which causes segmentation violation.
>>
>> Best,
>> Ehsan
>>
>>
>> On Sunday, July 3, 2016 at 4:32:16 PM UTC-5, bangerth wrote:
>>>
>>> On 07/03/2016 03:50 PM, Ehsan Esfahani wrote: 
>>> > Dear All, 
>>> > 
>>> > Greetings. I changed step-25 (minor changes) in order to solve my 
>>> problem. Now 
>>> > I want to change this code for parallel computation based on the 
>>> method 
>>> > mentioned in step-40. I got several errors and solved them one by one 
>>> but the 
>>> > following error: 
>>> > 
>>> > /Number of active cells: 1024 
>>> > //   Total number of cells: 1365 
>>> > //{[0,4224]}// 
>>> > //Time step #1; advancing to t = 0.1. 
>>>  > [...] 
>>> > //[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation 
>>> Violation, 
>>> > probably memory access out of range 
>>>  > [...] 
>>> > 
>>> > / 
>>> > / 
>>> > Eclipse gives me a backtrace in the following line of the code: 
>>> > /solver.solve (system_matrix, 
>>> completely_distributed_solution_update, 
>>> > system_rhs,/ 
>>> > /  preconditioner);/ 
>>> > I have no idea why I got this error. The code is running properly for 
>>> /fe(1) 
>>> > /and /n_global_refinements (4)/ but when I change them to /fe(2)/ and 
>>> > n_global_refinments (4) I got that error related to /Segmentation 
>>> Violation. 
>>> > /Do you know what's going on? Also, I have attached the code here . 
>>> Thanks in 
>>> > advance for your help. 
>>>
>>> Ehsan, 
>>>
>>> segmentation violations (SEGV) are typically easy to debug because you 
>>> can get 
>>> a backtrave in the debugger of the exact place where it happens, and you 
>>> can 
>>> then look at the local variables to see why this may have happened. Have 
>>> you 
>>> tried to run the program in a debugger and see what is going on? 
>>>
>>> Best 
>>>   W. 
>>>
>>> -- 
>>>  
>>> Wolfgang Bangerth   email:bang...@math.tamu.edu 
>>>  www: 
>>> http://www.math.tamu.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] External point forces on elastic bodies

[Hamed Babaei]

Dear Friends,

I want to solve a very simple elastic problem to obtain stress and 
displacement field, applying a point force on one of the external surfaces 
of a rectangular 2D domain.
Although I found some sort of external forces only in step 8, it is not a 
point force. I was wondering how to apply point forces in dealii.

Thanks,
Hamed

-- 
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: Thermoelastic Problem

[Hamed Babaei]

Hello Daniel,

Since although we can solve two elastic and heat equations separately, we 
still need to enter stress computed from elastic equation into the heat 
equation in every time step and keep updating stress field.
I was wondering when we have two different DoFHandlers for each of the 
equations, how to use, for instance, stress field stored on a quadrature 
point related to elastic problem DofHandler, in order to solve 
stress-induced heat equation problem. In other words, I don't know how to 
exchange information between two equation.

Thanks,
Hamed

On Friday, June 3, 2016 at 3:40:28 AM UTC-5, Daniel Arndt wrote:
>
> You're right, Hamed. There is FEValuesViews::Vector< dim, spacedim 
> >::symmetric_gradient [1] 
> and FEValuesViews::Vector< dim, spacedim >::gradient [2] instead. The 
> first should be exactly what you are looking for.
>
> Code that I have written normally uses multiple DoFHandlers if I solve for 
> multiple equations separately.
>
> Best,
> Daniel
>
> [1] 
> https://dealii.org/8.4.1/doxygen/deal.II/classFEValuesViews_1_1Vector.html#a4e5dfbb49d284a368acac6ef5dd4b2ef
>  
> <https://dealii.org/8.4.1/doxygen/deal.II/classFEValuesViews_1_1Vector.html>
> [2] 
> https://dealii.org/8.4.1/doxygen/deal.II/classFEValuesViews_1_1Vector.html#ad7b4df64147989229f566eff541ddebc
>  
> <https://dealii.org/8.4.1/doxygen/deal.II/classFEValuesViews_1_1Vector.html>
>
>
>
>
>
>
> Am Freitag, 3. Juni 2016 01:33:34 UTC+2 schrieb Hamed Babaei:
> Daniel,
>
> Since in the FEValuesViews::Vector< dim, spacedim > Class there isn't any 
> shape_grad_component function, I am not sure that I can update the 
> get_strain as I mentioned.
> I was wondering if you have already written a code for a thermoelastic 
> problem in dealii, solving heat equation and elastic equation separately.
>
> Best,
> Hamed
>
> On Thursday, June 2, 2016 at 3:31:48 PM UTC-5, Daniel Arndt wrote:
> Hamed,
>
>
> There shouldn't be a problem to do what you following step-22. In 
> particular, you can modify step-18::get_strain as you proposed.
> Did you face any difficulties so far?
>
>
> 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.