Re: [deal.II] Re: Issue regarding higher order element

2019-07-16 Thread Bruno Turcksin
Krishanu

Le mar. 16 juil. 2019 à 16:00, Krishanu Sen  a écrit :
>
> I am using deal.ii for a nonlinear elasticity problem. And the change in the 
> code that is causing this error message is:
>
>
>  template 
>   ElasticProblem::ElasticProblem ()
> :
> dof_handler (triangulation),
> fe (FE_Q(3), dim)
>   {}
No it's not where the error comes from. It is from somewhere in
Step8::ElasticProblem<1>::assemble_system_matrix() You need to find
which line produces the error.

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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/CAGVt9eOaE9g2sehRfaT-9B8JVq9pE79jdOV5uqMBF3LFzCQuaw%40mail.gmail.com.
For more options, visit https://groups.google.com/d/optout.


[deal.II] Re: Issue regarding higher order element

2019-07-16 Thread Krishanu Sen
Hello Bruno,

Here is the error message that I am getting:


An error occurred in line <388> of file 

 
in function
dealii::Tensor::value_type& dealii::Tensor::operator[](unsigned int) [with int rank_ = 2; int dim = 3; 
Number = double; dealii::Tensor::value_type = 
dealii::Tensor<1, 3, double>]
The violated condition was: 
i::assemble_system_matrix()
#1  ./nodal_growth_dyn_full: Step8::ElasticProblem<1>::run()
#2  ./nodal_growth_dyn_full: main


CMakeFiles/run.dir/build.make:57: recipe for target 'CMakeFiles/run' failed
make[3]: *** [CMakeFiles/run] Aborted
CMakeFiles/Makefile2:232: recipe for target 'CMakeFiles/run.dir/all' failed
make[2]: *** [CMakeFiles/run.dir/all] Error 2
CMakeFiles/Makefile2:239: recipe for target 'CMakeFiles/run.dir/rule' failed
make[1]: *** [CMakeFiles/run.dir/rule] Error 2
Makefile:183: recipe for target 'run' failed
make: *** [run] Error 2


I am using deal.ii for a nonlinear elasticity problem. And the change in 
the code that is causing this error message is:


 template 
  ElasticProblem::ElasticProblem ()
:
dof_handler (triangulation),
fe (FE_Q(3), dim)
  {}


The code runs fine when I am using FE_Q(2) or FE_Q(1) instead.

Thanks,
Krishanu

On Tuesday, 16 July 2019 14:22:13 UTC-5, Bruno Turcksin wrote:
>
> Krishanu,
>
> It's impossible to help you with so little information. Somewhere you set 
> something to the wrong size. Can you give us a backtrace? We need to see 
> the error message and the code that produces it.
>
> Best,
>
> Bruno
>
> On Tuesday, July 16, 2019 at 3:08:42 PM UTC-4, Krishanu Sen wrote:
>>
>> I am trying to run a code using higher order elements. The code runs fine 
>> when I am using linear [FE_Q(1)] or quadratic [FE_Q(2)] elements. 
>> But when I tried to use cubic elements [FE_Q(3)], the code stops with 
>> an error "ExcIndexRange", which says that the code is trying to access an 
>> index of a tensor that exceeds the limit of the tensor indices. I am not 
>> sure why would that happen just because of changing the polynomial degree 
>> of the element. I would appreciate any help or insight regarding the issue.
>>
>> Thanks,
>> Krishanu
>>
>

-- 
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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/0884e920-3ac7-4d63-84b5-d9993223e568%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.


[deal.II] Re: Issue regarding higher order element

2019-07-16 Thread Bruno Turcksin
Krishanu,

It's impossible to help you with so little information. Somewhere you set 
something to the wrong size. Can you give us a backtrace? We need to see 
the error message and the code that produces it.

Best,

Bruno

On Tuesday, July 16, 2019 at 3:08:42 PM UTC-4, Krishanu Sen wrote:
>
> I am trying to run a code using higher order elements. The code runs fine 
> when I am using linear [FE_Q(1)] or quadratic [FE_Q(2)] elements. 
> But when I tried to use cubic elements [FE_Q(3)], the code stops with 
> an error "ExcIndexRange", which says that the code is trying to access an 
> index of a tensor that exceeds the limit of the tensor indices. I am not 
> sure why would that happen just because of changing the polynomial degree 
> of the element. I would appreciate any help or insight regarding the issue.
>
> Thanks,
> Krishanu
>

-- 
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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/4f277d7a-396b-43f7-b887-32c850ae6944%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.


Re: [deal.II] Re: Compatibility of Petsc with step 18

2019-07-16 Thread Bruno Turcksin
Ramprasad,

On Tuesday, July 16, 2019 at 11:51:00 AM UTC-4, Ramprasad R wrote:
>
> The version of deal.ii is 9.0.1 and the version of PETSc is 3.10.5. I also 
> changed the config file of deal.ii to have MPI PETSc on.
>
What do you mean by "I changed the config file of deal.II"? You should not 
change any file generated by deal.II. When you compiled deal.II, did you 
use the option -DDEAL_II_WITH_PETSC=ON? If not, deal.II was not compiled 
with PETSc supports.

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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/0372d95b-ae52-457e-9400-9e8159469a62%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.


[deal.II] Issue regarding higher order element

2019-07-16 Thread Krishanu Sen
I am trying to run a code using higher order elements. The code runs fine 
when I am using linear [FE_Q(1)] or quadratic [FE_Q(2)] elements. 
But when I tried to use cubic elements [FE_Q(3)], the code stops with 
an error "ExcIndexRange", which says that the code is trying to access an 
index of a tensor that exceeds the limit of the tensor indices. I am not 
sure why would that happen just because of changing the polynomial degree 
of the element. I would appreciate any help or insight regarding the issue.

Thanks,
Krishanu

-- 
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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/7e068c9a-8a15-48b1-9f1e-2ed2453f33e8%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.


Re: [deal.II] Create a LinearOperator with a function instead of a matrix

2019-07-16 Thread Matthias Maier
A short comment:

On Tue, Jul 16, 2019, at 11:18 CDT, Matthias Maier  wrote:

> struct LeftVector {
> };

> struct RightVector {
> };

These two classes are of course just decoration (showing the minimal
interface a vector has to possess). There is usually no need to define
custom Vector clases.

-- 
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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/87d0iaar22.fsf%4043-1.org.
For more options, visit https://groups.google.com/d/optout.


Re: [deal.II] Create a LinearOperator with a function instead of a matrix

2019-07-16 Thread Matthias Maier
Hi,

the following example (which is a simplified version of test
lac/linear_operator_01) should get you started. In short: Simply create
an empty LinearOperator object (with the appropriate template
parameters) and populate the corresponding std::function objects.

Best,
Matthias



#include 
#include 

using namespace dealii;

struct LeftVector {
  typedef double value_type;
  value_type value;

  LeftVector =(value_type new_value)
  {
value = new_value;
return *this;
  }
  LeftVector *=(value_type scale)
  {
value *= scale;
return *this;
  }
  LeftVector /=(value_type scale)
  {
value /= scale;
return *this;
  }
  LeftVector +=(const LeftVector )
  {
value += u.value;
return *this;
  }
  int size() const
  {
return 1;
  }
  std::size_t memory_consumption() const
  {
return 1;
  }
};

struct RightVector {
  typedef double value_type;
  value_type value;

  RightVector =(value_type new_value)
  {
value = new_value;
return *this;
  }
  RightVector *=(value_type scale)
  {
value *= scale;
return *this;
  }
  RightVector /=(value_type scale)
  {
value /= scale;
return *this;
  }
  RightVector +=(const RightVector )
  {
value += u.value;
return *this;
  }
  int size() const
  {
return 1;
  }
  std::size_t memory_consumption() const
  {
return 1;
  }
};

int main()
{
  LinearOperator multiply2;

  multiply2.vmult = [](LeftVector , const RightVector ) {
v.value = 2 * u.value;
  };
  multiply2.vmult_add = [](LeftVector , const RightVector ) {
v.value += 2 * u.value;
  };
  multiply2.Tvmult = [](RightVector , const LeftVector ) {
v.value = 2 * u.value;
  };
  multiply2.Tvmult_add = [](RightVector , const LeftVector ) {
v.value += 2 * u.value;
  };
  multiply2.reinit_range_vector = [](LeftVector &, bool omit_zeroing_values) {
// do nothing
  };
  multiply2.reinit_domain_vector = [](RightVector &, bool omit_zeroing_values) {
// do nothing
  };

  // Small test:

  RightVector u = {4.};
  LeftVector v = {0.};

  multiply2.vmult(v, u);
  std::cout << "2 * " << u.value << " = " << v.value << std::endl;

  multiply2.vmult_add(v, u);
  std::cout << "... + 2 * " << u.value << " = " << v.value << std::endl;

  return 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 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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/87h87mar6y.fsf%4043-1.org.
For more options, visit https://groups.google.com/d/optout.


[deal.II] Create a LinearOperator with a function instead of a matrix

2019-07-16 Thread 'Maxi Miller' via deal.II User Group
I tried to follow the examples for setting up a LinearOperator (such as 
shown in Step-22 and Step-20), but instead of a matrix I wanted to provide 
my own function which should serve as vmult-function (i.e. a matrix-free 
linear operator). Is that possible? My initial test was the following code:

  class approximator : public Subscriptor 
{ 
public: 
   approximator(std::function&, 
Vector&)> 
residual_function, 
const Vector _point) : 
residual_function(residual_function)
{ 
u = evaluation_point; 
u.add(1); 
epsilon_val = u.l1_norm() * 1e-6; 
residual_function(u, residual_0); 
}; 

void vmult(Vector , const Vector ) const; 

private: 
std::function&, 
   Vector&)> residual_function; 
double epsilon_val; 
Vector u, residual_0, residual_1; 
};


const auto op_M = linear_operator(jacobian_approximation(std::bind(&
MinimalSurfaceProblem::compute_residual,
   this,
   std::
placeholders::_1,
   std::
placeholders::_2), 

 present_solution));

but that resulted in
error: use of deleted function ‘dealii::LinearOperator dealii::linear_operator(Matrix&&) [with Range = dealii::Vector
; Domain = dealii::Vector; Payload = dealii::internal::
LinearOperatorImplementation::EmptyPayload; Matrix = Step15::
jacobian_approximation;  = void]’
  389 | 
 present_solution));
Are there other ways?
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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/1cf547d8-4d9b-4826-b0c8-85c2bd2d2f3f%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.


Re: [deal.II] Re: Compatibility of Petsc with step 18

2019-07-16 Thread Ramprasad R
Dear Daniel,
The version of deal.ii is 9.0.1 and the version of PETSc is 3.10.5. I also
changed the config file of deal.ii to have MPI PETSc on.

Thank you.

Best regards,
Ramprasad

On Tue, Jul 16, 2019 at 5:12 PM Daniel Arndt  wrote:

> Ramprasad,
>
> Let's first rule out some simple problems. Which version of deal.II are
> you using? Did you compile deal.II with MPI and PETSc support?
> We only define the namecpacePETScWrappers when DEAL_II_WITH_PETSC=ON. Can
> you provide us with the content of the detailed.log file in your build
> directory?
>
> 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.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/1b94eabe-7cd4-4b79-8ffb-207234fb026c%40googlegroups.com
> 
> .
> For more options, visit https://groups.google.com/d/optout.
>

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


Re: [deal.II] Implicit stepping methods from the TimeStepping-namespace for large sparse matrices

2019-07-16 Thread 'Maxi Miller' via deal.II User Group
If I understand it correctly, that is equivalent to solving Ax = b, with A 
and b given, and x returned?

Am Montag, 15. Juli 2019 22:03:46 UTC+2 schrieb Daniel Arndt:
>
> Maxi,
>
> The interface for TimeStepping::evolve_one_time_step is:
>
> virtual double
>   evolve_one_time_step(
>   std::vector>
>   & F,
>   std::vector   VectorType(const double, const double, const VectorType &)>> _inverse,
>   double t,
>   double delta_t,
>   VectorType & y);
>
> and does not depend on UMFPACK at all. You just need to provide a 
> function-type object that can be used for evaluating the inverse of the 
> Jacobians and the right-hand side.
> You might also want to have a look at the test base/time_stepping_01 (
> https://github.com/dealii/dealii/blob/master/tests/base/time_stepping_01.cc
> ).
>
> Best,
> Daniel
>
>
> Am Mo., 15. Juli 2019 um 15:04 Uhr schrieb 'Maxi Miller' via deal.II User 
> Group >:
>
>> As far as I understand all implicit time-stepping methods in the 
>> TimeStepping-namespace take a function which invert the jacobian matrix and 
>> the mass matrix (combined). This is done using a sparse solver (UMFPACK). 
>> Is it that still possible for a larger system (10kk DoFs), or do I have to 
>> resort to other methods, after inverting the sparse matrix will result in a 
>> dense matrix (which likely will overflow my memory)?
>>
>> -- 
>> 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 dea...@googlegroups.com .
>> To view this discussion on the web visit 
>> https://groups.google.com/d/msgid/dealii/ccfe4a75-775d-471b-b69e-4a886c2e443d%40googlegroups.com
>>  
>> 
>> .
>> For more options, visit https://groups.google.com/d/optout.
>>
>

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


[deal.II] Re: Compatibility of Petsc with step 18

2019-07-16 Thread Daniel Arndt
Ramprasad,

Let's first rule out some simple problems. Which version of deal.II are you 
using? Did you compile deal.II with MPI and PETSc support?
We only define the namecpacePETScWrappers when DEAL_II_WITH_PETSC=ON. Can 
you provide us with the content of the detailed.log file in your build 
directory?

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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/1b94eabe-7cd4-4b79-8ffb-207234fb026c%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.


[deal.II] Compatibility of Petsc with step 18

2019-07-16 Thread Ramprasad R
Hello All,

I am getting an error when I try to compile the step-18. It states 
"PETScWrappers does not name a type". I wanted to know if anyone else has 
faced the same problem. And if there is a solution to this.

Thank you. 

-- 
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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/2d25a910-8d32-4268-a8aa-818ff5416670%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.


[deal.II] Re: Porting tutorials to PETSc from Trilinos

2019-07-16 Thread Bruno Turcksin
Franco,

On Tuesday, July 16, 2019 at 3:54:34 AM UTC-4, Franco Milicchio wrote:
>
>
> Can I ask why there are some discrepancies between those two wrappers? 
> Apart from an obvious one, that they're different libraries, but uniformity 
> may help a lot IMO.
>
Because the wrappers have been written over time by different people who 
were not aware of how exactly the other wrapper works (most people use 
either Trilinos or PETSc not both). We are trying to get the interfaces to 
be more similar but it is very difficult because any change in the 
interface will break many codes.

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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/533ea4e2-44fd-43f0-8bd7-93b77bdbd6ed%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.


[deal.II] Re: Porting tutorials to PETSc from Trilinos

2019-07-16 Thread Franco Milicchio

>
>
> step-40 tries to allow switching between Trilinos and PETSc as easy as 
> possible by using the deal.II/lac/generic_linear_algebra.h header and the 
> defining the LA::MPI namespace.
>

Thanks for your answer, Daniel. Even with that example I find it a little 
harder than I expected to use classes from PETSc instead of Trilinos. For 
instance:

const LinearSolvers::InverseMatrix
mp_inverse(stokes_preconditioner_matrix.block(1, 1), 
*Mp_preconditioner);

Cannot compile, I am trying to use a 
static_cast(stokes_preconditioner_matrix.block(1, 
1)) now, but I don't think it will be a good idea.
 
Another weak point in my porting is

template 
BlockSchurPreconditioner::
  BlockSchurPreconditioner(
const PETScWrappers::MPI::BlockSparseMatrix ,
const InverseMatrix ,
const PreconditionerTypeA &Apreconditioner)
  : stokes_matrix()
  , m_inverse()
  , a_preconditioner(Apreconditioner)
  , tmp(complete_index_set(stokes_matrix->block(1, 1).m()))
{}

where tmp (a mutable PETScWrappers::MPI::Vector) cannot be initialized. 
Adding MPI_COMM_WORLD make it compile, though.

Again with matrices, I have made the following changes:

//temperature_matrix.copy_from(temperature_mass_matrix);
temperature_matrix = 0;
temperature_matrix.add(1.0, temperature_mass_matrix);

I hope the change works...

Can I ask why there are some discrepancies between those two wrappers? 
Apart from an obvious one, that they're different libraries, but uniformity 
may help a lot IMO.

 
>
>> [...]
>> std::shared_ptr Amg_preconditioner;
>>
>  
> This looks inappropriate. Probably it should rather be 
> PETScWrappers::PreconditionBoomerAMG;
>


Yes, I know :) I am first trying to make it compile, then I will try to use 
a good preconditioner (Boomer requires Hypre, and it needs MPI, so right 
now it is not an option).

Thanks!
Franco

-- 
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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/0bce36bd-19a1-4f62-8be8-75eb3965e7ce%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.