Re: [deal.II] Linear operators - TrilinosWrappers

[Jean-Paul Pelteret]

Dear Felipe,

It might be that you need to set up the preconditioner operator with an
exemplar matrix (since the identity operator doesn't know the size of the
range and domain that it's working on).

If that's not the issue then could you please try to reproduce this as a
minimal example and share it with us? The program doesn't need to produce
any meaningful result, but it would be good if it shows both the working
scenario and the problematic one. That way we could investigate further and
try to figure out what's going wrong here.

Best,
Jean-Paul

Sent from my mobile device. Please excuse my brevity and any typos.

On Wed, 16 Jun 2021, 08:50 Juan Felipe Giraldo, 
wrote:

> Dear Jean-Paul,
>
> Thank you for your reply and the complete information you provide me.
> Effectively, declaring the preconditioner (out-of-line) was one problem,
> but the inverse_operator function is still not matching.
>
> What I just realized is that when I compute the operator as a
> multiplication of a linear operator, an inverse operator and a linear
> operator ( for instance, when I compute the Schur complement )
>
> const auto op_S = op_BT * op_M_inv * op_B;
>
> this operator becomes to the type
> 'dealii::TrilinosWrappers::internal::LinearOperatorImplementation::TrilinosPayload'
> ,
>
> but when I compute the operator as a multiplication of other linear
> operators directly (for instance, when I compute the Schur preconditioner,
> following the procedure in step 20 )
>
> const auto op_pre = linear_operator(prec_M);
> const auto op_aS = op_BT * op_pre * op_B;
>
> this operator becomes to the type
> 'dealii::internal::LinearOperatorImplementation::EmptyPayload&'.
>
> So, If I use the inverse_operator function with the first operator (op_S),
> I can obtain the inverse without any problem,
> (const auto op_S_inv = inverse_operator(op_S, solver_S, preconditioner_S );
>
> But if I do it with the second case, it doesn't work because the functions
> are not matching.
> (const auto op_S_inv = inverse_operator(opa_S, solver_S, preconditioner_S
> );
>
> I am not sure what is happening because all the operators are declared as
> "LA::MPI::Vector".
> If you have any suggestion, would be greatly appreciated.
>
> Thank you so much,
>
> Regards,
> Felipe Giraldo
>
>
>
> El martes, 15 de junio de 2021 a las 19:16:54 UTC+8, Jean-Paul Pelteret
> escribió:
>
>> Hi again Feilpe,
>>
>> Regarding the lack of documentation, I’ve opened an issue on Github to
>> track this. You can find that here:
>> https://github.com/dealii/dealii/issues/12466
>>
>> Best,
>> Jean-Paul
>>
>> On 15. Jun 2021, at 12:57, Jean-Paul Pelteret  wrote:
>>
>> Hi Feilpe,
>>
>> Firstly, I agree that the documentation is very light on details on how
>> to use the linear operators with Trilinos linear algebra types. We can
>> definitely improve this, and any contributions to that effect would be
>> greatly appreciated!
>>
>> Right now, I can direct you to a few of the tests that use Trilinos LA in
>> conjunction with the inverse operator, so that you can compare what you
>> have and figure out what the problematic differences are. There is this
>> one, for instance
>>
>> https://github.com/dealii/dealii/blob/master/tests/lac/linear_operator_12a.cc#L323-L341
>> that looks like a similar setup to yours, and
>>
>> https://github.com/dealii/dealii/blob/master/tests/lac/schur_complement_04.cc#L126-L137
>> that uses a Trilinos::SolverCG (as opposed to deal.II’s solver).
>>
>> Comparing to both of these, I think that the important point might be
>> that the preconditioner must be declared (out-of-line) before the
>> inverse_operation() function is called. The linear operators typically
>> expect the lifetime of the LA objects to exceed that of the operators, and
>> so you have to first create the matrix or preconditioner and then pass it
>> to the linear operators. This is similar to what you’d done when setting up 
>> op_M,
>> for example. The case of passing in a deal::PreconditionerIdentity() for
>> serial operations is a special case, and I don’t think that we’ve
>> duplicated that for TrilinosWrappers:: PreconditionerIdentity(). Maybe
>> that could be improved too.
>>
>> I hope that with this single change you’d be able to get your
>> program running. If not, then please do let us know so that we can try to
>> help further.
>>
>> Best,
>> Jean-Paul
>>
>>
>> On 14. Jun 2021, at 19:18, Juan Felipe Giraldo 
>> wrote:
>>
>> Hello everyone,
>>
>> I have implemented a residual-minimization framework that somehow is
>> similar to DPG. I want to extend my results by using parallelization using
>> MPI with PETSc or Trilinos.
>> So far, I have solved the saddle point problem using the Schur complement
>> exactly how it is described in step 20. Now, I am trying to replicate
>> exactly the same solver but using the MPI wrappers and the linear operators.
>>
>> The problem is that when I am trying to implement the inverse_operator to
>> compute the Preconditioner of the Schur complement, I g

Re: [deal.II] Linear operators - TrilinosWrappers

[Juan Felipe Giraldo]

Dear Jean-Paul,

Thank you for your reply,

I already set up the Identity preconditioner with a block matrix, but the 
problem remains the same. I am not sure this is the problem because I 
didn't need to set up the identity operator to any matrix in the scenario 
where I could invert. I think my problem is in the Schur preconditioner 
operator before to invert.
I attach the problematic code and the "working code" I am coding. The only 
difference between them is the line 1251:

(const auto op_S_inv = inverse_operator(opa_S, solver_S, preconditioner_S 
); //(problematic case)
(const auto op_S_inv = inverse_operator(op_S, solver_S, preconditioner_S 
);  //(Case I could invert)

(Please excuse me for the extension of the code, but the solver function 
where I have the issue is too short. In the case you consider it necessary, 
I could code a shorter one)

Thank you so much,

Regards,
Felipe



El miércoles, 16 de junio de 2021 a las 15:45:27 UTC+8, Jean-Paul Pelteret 
escribió:

> Dear Felipe,
>
> It might be that you need to set up the preconditioner operator with an 
> exemplar matrix (since the identity operator doesn't know the size of the 
> range and domain that it's working on).
>
> If that's not the issue then could you please try to reproduce this as a 
> minimal example and share it with us? The program doesn't need to produce 
> any meaningful result, but it would be good if it shows both the working 
> scenario and the problematic one. That way we could investigate further and 
> try to figure out what's going wrong here. 
>
> Best,
> Jean-Paul
>
> Sent from my mobile device. Please excuse my brevity and any typos.
>
> On Wed, 16 Jun 2021, 08:50 Juan Felipe Giraldo,  
> wrote:
>
>> Dear Jean-Paul,
>>
>> Thank you for your reply and the complete information you provide me. 
>> Effectively, declaring the preconditioner (out-of-line) was one problem, 
>> but the inverse_operator function is still not matching. 
>>
>> What I just realized is that when I compute the operator as a 
>> multiplication of a linear operator, an inverse operator and a linear 
>> operator ( for instance, when I compute the Schur complement )
>>
>> const auto op_S = op_BT * op_M_inv * op_B;
>>
>> this operator becomes to the type 
>> 'dealii::TrilinosWrappers::internal::LinearOperatorImplementation::TrilinosPayload'
>>  
>> ,
>>
>> but when I compute the operator as a multiplication of other linear 
>> operators directly (for instance, when I compute the Schur preconditioner, 
>> following the procedure in step 20 )
>>
>> const auto op_pre = linear_operator(prec_M);
>> const auto op_aS = op_BT * op_pre * op_B;
>>
>> this operator becomes to the type 
>> 'dealii::internal::LinearOperatorImplementation::EmptyPayload&'.
>>
>> So, If I use the inverse_operator function with the first operator 
>> (op_S), I can obtain the inverse without any problem,
>> (const auto op_S_inv = inverse_operator(op_S, solver_S, preconditioner_S 
>> );
>>
>> But if I do it with the second case, it doesn't work because the 
>> functions are not matching. 
>> (const auto op_S_inv = inverse_operator(opa_S, solver_S, preconditioner_S 
>> );
>>
>> I am not sure what is happening because all the operators are declared 
>> as  "LA::MPI::Vector".
>> If you have any suggestion, would be greatly appreciated. 
>>
>> Thank you so much,
>>
>> Regards,
>> Felipe Giraldo
>>
>>
>>
>> El martes, 15 de junio de 2021 a las 19:16:54 UTC+8, Jean-Paul Pelteret 
>> escribió:
>>
>>> Hi again Feilpe,
>>>
>>> Regarding the lack of documentation, I’ve opened an issue on Github to 
>>> track this. You can find that here: 
>>> https://github.com/dealii/dealii/issues/12466
>>>
>>> Best,
>>> Jean-Paul
>>>
>>> On 15. Jun 2021, at 12:57, Jean-Paul Pelteret  
>>> wrote:
>>>
>>> Hi Feilpe,
>>>
>>> Firstly, I agree that the documentation is very light on details on how 
>>> to use the linear operators with Trilinos linear algebra types. We can 
>>> definitely improve this, and any contributions to that effect would be 
>>> greatly appreciated!
>>>
>>> Right now, I can direct you to a few of the tests that use Trilinos LA 
>>> in conjunction with the inverse operator, so that you can compare what you 
>>> have and figure out what the problematic differences are. There is this 
>>> one, for instance
>>>
>>> https://github.com/dealii/dealii/blob/master/tests/lac/linear_operator_12a.cc#L323-L341
>>> that looks like a similar setup to yours, and
>>>
>>> https://github.com/dealii/dealii/blob/master/tests/lac/schur_complement_04.cc#L126-L137
>>> that uses a Trilinos::SolverCG (as opposed to deal.II’s solver).
>>>
>>> Comparing to both of these, I think that the important point might be 
>>> that the preconditioner must be declared (out-of-line) before the 
>>> inverse_operation() function is called. The linear operators typically 
>>> expect the lifetime of the LA objects to exceed that of the operators, and 
>>> so you have to first create the matrix or preconditioner and then pass it 
>>> to

[deal.II] Neumann data and Kelly estimator: can't pass std::map argument

[bob bill]

Hi everyone,

I'm trying to change a bit the step-7 (Helmotz Problem 

)

In particular, I'm going to change the estimate function
KellyErrorEstimator::estimate 

(
dof_handler,
QGauss 
(fe->degree 
+ 1),
std::maphttps://www.dealii.org/current/doxygen/deal.II/classunsigned_01int.html>, 
const Function 
 *>(),
solution,
estimated_error_per_cell);

since I want to encode the fact that we have Neumann data. To do so, I 
created a function that describes my Neumann data 


template
class NeumannData: public Function { 
virtual double value(const Point &p, const unsigned int) const 
override; 
}; 

template 
double NeumannData::value(const Point &p, const unsigned int) 
const {return p[0] + p[1];}


and hence my refine function becomes (in red what I changed w.r.t the 
original tutorial program)

template

void step7::refine_grid() {
const NeumannData neumann_func;

switch (refinement_mode) {
case global_refinement: {
triangulation.refine_global(1);

break;

}

case adaptive_refinement: {

Vector estimated_error_per_cell(triangulation.n_active_cells());

KellyErrorEstimator::estimate(dof_handler,

QGauss(fe->degree + 1),

std::map *>(1,&neumann_func),

estimated_error_per_cell);
GridRefinement::refine_and_coarsen_fixed_number(triangulation,

estimated_error_per_cell, 0.3, 0.03);

triangulation.execute_coarsening_and_refinement();
break;

}


default: {

Assert(false, ExcNotImplemented());

}

}

}

*The boundary indicator for Neumann data is 1, and I gave a pointer to a 
constant function as argument.  **Unfortunately, the compiler says: "*no 
matching constructor for initialization of 'std::map *>' (aka 'map *>') 
std::map *>(1,&neumann_func),*"*

*What am I missing?*


Thanks,

Bob

-- 
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/53039ccf-c4f3-4b3b-b6f7-e8f3590b1fcfn%40googlegroups.com.

Re: [deal.II] Linear operators - TrilinosWrappers

[Jean-Paul Pelteret]

Dear Felipe,

Changing this
const auto op_pre = linear_operator(prec_M);
to this
const auto op_pre = linear_operator(M, prec_M);
allowed your program "FEMDG-Dar-par_tofix" to compile for me. If you have a 
look at this documentation for this variant of linear_operator
https://dealii.org/current/doxygen/deal.II/group__LAOperators.html#ga14dbc8c2c27ea3fd45576528a891c6e2
 

you'll see that it specifically mentions its use to encapsulate preconditioners 
— this is what I was referring to in my previous message. It’s not a 
particularly easy problem to diagnose, so I’ll give some thought on what (if 
anything) we can do to programatically improve this situation on our side. In 
the mean time, I've added a note to the GitHub issue to detail this further.

(BTW. You also probably want to change this
const auto op_M  = linear_operator(M, prec_M);
to this
const auto op_M = linear_operator(system_matrix.block(0, 0));
if you are to use it, as op_M then acts like a preconditioner for M, which I 
don’t think is what you want.)

I didn’t try to run the code after checking that it compiled, but I hope that 
this now sorts things out for you. If not, just let us know.

Best,
Jean-Paul

> On 16. Jun 2021, at 12:19, Juan Felipe Giraldo  wrote:
> 
> Dear Jean-Paul,
> 
> Thank you for your reply,
> 
> I already set up the Identity preconditioner with a block matrix, but the 
> problem remains the same. I am not sure this is the problem because I didn't 
> need to set up the identity operator to any matrix in the scenario where I 
> could invert. I think my problem is in the Schur preconditioner operator 
> before to invert.
> I attach the problematic code and the "working code" I am coding. The only 
> difference between them is the line 1251:
> 
> (const auto op_S_inv = inverse_operator(opa_S, solver_S, preconditioner_S ); 
> //(problematic case)
> (const auto op_S_inv = inverse_operator(op_S, solver_S, preconditioner_S );  
> //(Case I could invert)
> 
> (Please excuse me for the extension of the code, but the solver function 
> where I have the issue is too short. In the case you consider it necessary, I 
> could code a shorter one)
> 
> Thank you so much,
> 
> Regards,
> Felipe
> 
> 
> 
> El miércoles, 16 de junio de 2021 a las 15:45:27 UTC+8, Jean-Paul Pelteret 
> escribió:
> Dear Felipe,
> 
> It might be that you need to set up the preconditioner operator with an 
> exemplar matrix (since the identity operator doesn't know the size of the 
> range and domain that it's working on).
> 
> If that's not the issue then could you please try to reproduce this as a 
> minimal example and share it with us? The program doesn't need to produce any 
> meaningful result, but it would be good if it shows both the working scenario 
> and the problematic one. That way we could investigate further and try to 
> figure out what's going wrong here. 
> 
> Best,
> Jean-Paul
> 
> Sent from my mobile device. Please excuse my brevity and any typos.
> 
> On Wed, 16 Jun 2021, 08:50 Juan Felipe Giraldo,  > wrote:
> Dear Jean-Paul,
> 
> Thank you for your reply and the complete information you provide me. 
> Effectively, declaring the preconditioner (out-of-line) was one problem, but 
> the inverse_operator function is still not matching. 
> 
> What I just realized is that when I compute the operator as a multiplication 
> of a linear operator, an inverse operator and a linear operator ( for 
> instance, when I compute the Schur complement )
> 
> const auto op_S = op_BT * op_M_inv * op_B;
> 
> this operator becomes to the type 
> 'dealii::TrilinosWrappers::internal::LinearOperatorImplementation::TrilinosPayload'
>  ,
> 
> but when I compute the operator as a multiplication of other linear operators 
> directly (for instance, when I compute the Schur preconditioner, following 
> the procedure in step 20 )
> 
> const auto op_pre = linear_operator(prec_M);
> const auto op_aS = op_BT * op_pre * op_B;
> 
> this operator becomes to the type 
> 'dealii::internal::LinearOperatorImplementation::EmptyPayload&'.
> 
> So, If I use the inverse_operator function with the first operator (op_S), I 
> can obtain the inverse without any problem,
> (const auto op_S_inv = inverse_operator(op_S, solver_S, preconditioner_S );
> 
> But if I do it with the second case, it doesn't work because the functions 
> are not matching. 
> (const auto op_S_inv = inverse_operator(opa_S, solver_S, preconditioner_S );
> 
> I am not sure what is happening because all the operators are declared as  
> "LA::MPI::Vector".
> If you have any suggestion, would be greatly appreciated. 
> 
> Thank you so much,
> 
> Regards,
> Felipe Giraldo
> 
> 
> 
> El martes, 15 de junio de 2021 a las 19:16:54 UTC+8, Jean-Paul Pelteret 
> escribió:
> Hi again Feilpe,
> 
> Regarding the lack of documentation, I’ve opened an issue on Github to track 
> this. You can find that here: https://github.com/dealii

Re: [deal.II] Linear operators - TrilinosWrappers

[Juan Felipe Giraldo]

Dear Jean-Paul,

Finally it is working! Thank you for your help and clear explanation!   Now
it is a bit more clear for me how these linear operators work here :).

(Regarding the op_M operator,  I was not using that variable because I
declared M directly as the linear operator of the block(0,0). But thanks to
let me know to avoid problems!)

kind regards,
Felipe


On Wed, Jun 16, 2021 at 7:22 PM Jean-Paul Pelteret 
wrote:

> Dear Felipe,
>
> Changing this
> const auto op_pre = linear_operator(prec_M);
> to this
> const auto op_pre = linear_operator(M, prec_M);
> allowed your program "FEMDG-Dar-par_tofix" to compile for me. If you have
> a look at this documentation for this variant of linear_operator
>
> https://dealii.org/current/doxygen/deal.II/group__LAOperators.html#ga14dbc8c2c27ea3fd45576528a891c6e2
> you'll see that it specifically mentions its use to encapsulate
> preconditioners — this is what I was referring to in my previous message.
> It’s not a particularly easy problem to diagnose, so I’ll give some thought
> on what (if anything) we can do to programatically improve this situation
> on our side. In the mean time, I've added a note to the GitHub issue to
> detail this further.
>
> (BTW. You also probably want to change this
> const auto op_M = linear_operator(M, prec_M);
> to this
> const auto op_M = linear_operator(system_matrix.block(0,
> 0));
> if you are to use it, as op_M then acts like a preconditioner for M, which
> I don’t think is what you want.)
>
> I didn’t try to run the code after checking that it compiled, but I hope
> that this now sorts things out for you. If not, just let us know.
>
> Best,
> Jean-Paul
>
> On 16. Jun 2021, at 12:19, Juan Felipe Giraldo 
> wrote:
>
> Dear Jean-Paul,
>
> Thank you for your reply,
>
> I already set up the Identity preconditioner with a block matrix, but the
> problem remains the same. I am not sure this is the problem because I
> didn't need to set up the identity operator to any matrix in the scenario
> where I could invert. I think my problem is in the Schur preconditioner
> operator before to invert.
> I attach the problematic code and the "working code" I am coding. The only
> difference between them is the line 1251:
>
> (const auto op_S_inv = inverse_operator(opa_S, solver_S, preconditioner_S
> ); //(problematic case)
> (const auto op_S_inv = inverse_operator(op_S, solver_S, preconditioner_S
> );  //(Case I could invert)
>
> (Please excuse me for the extension of the code, but the solver function
> where I have the issue is too short. In the case you consider it necessary,
> I could code a shorter one)
>
> Thank you so much,
>
> Regards,
> Felipe
>
>
>
> El miércoles, 16 de junio de 2021 a las 15:45:27 UTC+8, Jean-Paul Pelteret
> escribió:
>
>> Dear Felipe,
>>
>> It might be that you need to set up the preconditioner operator with an
>> exemplar matrix (since the identity operator doesn't know the size of the
>> range and domain that it's working on).
>>
>> If that's not the issue then could you please try to reproduce this as a
>> minimal example and share it with us? The program doesn't need to produce
>> any meaningful result, but it would be good if it shows both the working
>> scenario and the problematic one. That way we could investigate further and
>> try to figure out what's going wrong here.
>>
>> Best,
>> Jean-Paul
>>
>> Sent from my mobile device. Please excuse my brevity and any typos.
>>
>> On Wed, 16 Jun 2021, 08:50 Juan Felipe Giraldo, 
>> wrote:
>>
>>> Dear Jean-Paul,
>>>
>>> Thank you for your reply and the complete information you provide me.
>>> Effectively, declaring the preconditioner (out-of-line) was one problem,
>>> but the inverse_operator function is still not matching.
>>>
>>> What I just realized is that when I compute the operator as a
>>> multiplication of a linear operator, an inverse operator and a linear
>>> operator ( for instance, when I compute the Schur complement )
>>>
>>> const auto op_S = op_BT * op_M_inv * op_B;
>>>
>>> this operator becomes to the type
>>> 'dealii::TrilinosWrappers::internal::LinearOperatorImplementation::TrilinosPayload'
>>> ,
>>>
>>> but when I compute the operator as a multiplication of other linear
>>> operators directly (for instance, when I compute the Schur preconditioner,
>>> following the procedure in step 20 )
>>>
>>> const auto op_pre = linear_operator(prec_M);
>>> const auto op_aS = op_BT * op_pre * op_B;
>>>
>>> this operator becomes to the type
>>> 'dealii::internal::LinearOperatorImplementation::EmptyPayload&'.
>>>
>>> So, If I use the inverse_operator function with the first operator
>>> (op_S), I can obtain the inverse without any problem,
>>> (const auto op_S_inv = inverse_operator(op_S, solver_S, preconditioner_S
>>> );
>>>
>>> But if I do it with the second case, it doesn't work because the
>>> functions are not matching.
>>> (const auto op_S_inv = inverse_operator(opa_S, solver_S,
>>> preconditioner_S );
>>>
>>> I am not sure what is happening because all the 

[deal.II] Re: online deal.II workshop: Friday, June 18

[Timo Heister]

Hi all,

This is a quick reminder that our one-day workshop will happen this Friday.
We hope to see you there and celebrate the 9.3 release with us!

The information on how to join the meeting via Zoom is now available at
https://dealii.org/workshop-2021/

Best,
Timo







On Wed, Jun 2, 2021 at 11:35 AM Timo Heister  wrote:
>
> Hi all,
>
> We would like to announce a one-day deal.II workshop on June 18, 2021
> with several talks about recent developments in the library, the 9.3
> release, and interesting user projects. The talks will be available
> live using zoom and available as a video at a later time.
>
> For more information including a schedule, please see:
> https://dealii.org/workshop-2021/
>
> Information on how to join will be made available closer to the date.
> We hope to see you there!
>
> Best,
> Timo and the deal.II developers
>
> --
> Timo Heister
> http://www.math.clemson.edu/~heister/



--
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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/CAMRj59GboArcKgkPTgPETRXcOerNDOgELt%2BvdNGSfEFLR%2B52%2Bg%40mail.gmail.com.

[deal.II] Local Raviart-Thomas space

[Giselle Sosa Jones]

Hello,

I am trying to implement a projection of a DG solution to a Raviart-Thomas 
space. For this, I need a local space Q_{p-1,p}(K) \times Q_{p,p-1}(K). I 
have been using FE_DGRaviartThomas of degree p-1, and then I 
use shape_value_component(i, q, 1) when I want to use the x component of 
the test function, and shape_value_component(i,q,0) when I want to use the 
y component of the test function. I have a bug in my code and I am 
wondering if it is coming from the way I am creating this local space. Does 
this construction of the space make sense? Is there maybe a better way of 
doing it?

Many thanks.

Giselle

-- 
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/5cb3de29-b3ba-4c81-9339-543253d67ee4n%40googlegroups.com.

Re: [deal.II] Re: Integrated material and spatial traction forces on boundary not equal

[Michael Lee]

Hi Alex,
I just got a chance to have Trilinos installed but the same error occurs 
regarding *EnergyFunctional*. I got the following notice when trying to 
test *step-31*:

Error! This tutorial requires a deal.II library that was configured with
  the following options:

  DEAL_II_WITH_TRILINOS = ON

  However, the deal.II library found at /usr/local was configured with these
  options

  DEAL_II_WITH_TRILINOS = OFF

  which conflict with the requirements.


-- Configuring incomplete, errors occurred!  

Is there a way to let deal.II know Trillinos is installed? Do I have to 
reinstall deal.ii?
Or, can I compile the code without using Trillinos with some modification?

Thanks!
Michael

On Friday, June 11, 2021 at 4:04:20 AM UTC-6 alexanderc...@gmail.com wrote:

> Hi Michael,
>
> I think that is the first point in my code that something from the 
> automatic differentiation module is referred to, which is a part of 
> trilinos. Was your version of deal.ii compiled with trilinos, and if so, 
> was your version of trilinos compiled with sacado?
>
> To answer the question I posed, I did some more tests and found that 
> refining in the width and height of the beam has no effect on the 
> convergence of the shear force. I went down to just a single division in z 
> and 2 in y (in other words there are just two faces when looking down the 
> end of the beam), and could get the same level of agreement with beam 
> theory with one order of magnitude fewer degrees of freedom. I guess I am 
> still mildly surprised how refined I had to made the x direction, but 
> perhaps if I also used adaptive refinement things would further improve.
>
> Best,
> Alex
>
> On Thursday, June 10, 2021 at 6:20:58 p.m. UTC+2 lian...@gmail.com wrote:
>
>> Alex, 
>>
>>  
>>
>> Thank you for sharing your codes. I have some compiling errors relating 
>> to “EnergyFunctional’: 
>>
>> solve_ring_nonlinear.cpp:520:43: error: ‘*EnergyFunctional’* in 
>> namespace ‘dealii::Differentiation::AD’ *does not name a template type*
>>
>>   520 | using ADHelper = Differentiation::AD::EnergyFunctional<
>>
>>   |   ^~~~
>>
>> Can you help me with that?
>>
>>  
>>
>> For the large error, I noticed you’re using linear element. I encountered 
>> the same large error when comparing it with Abaqus with FE_Q(1). But the 
>> error came down with less grids when I used higher finite element 
>> FE_Q(2). I remember the deflection of beam is a cubic function of 
>> coordinate. You may try that see if it improves.
>>
>>  
>>
>> Best,
>>
>> Michael
>>
>>  
>>
>> *From: *Alex Cumberworth
>> *Sent: *Wednesday, June 9, 2021 9:54 AM
>> *To: *deal.II User Group
>> *Subject: *[deal.II] Re: Integrated material and spatial traction forces 
>> on boundary not equal
>>
>>  
>>
>> Hello,
>>
>>  
>>
>> I have attached the most recent version of my code here. I have tried to 
>> make setting boundary conditions in the parameter file more convenient for 
>> myself; you can set boundary domains, boundary conditions that use these 
>> domains, and stages if the displacement is large. However, there are not 
>> many comments, so you may want to just remove this part for your own 
>> purposes.
>>
>>  
>>
>> However, I have been a bit surprised in my comparisons at how fine a mesh 
>> is required to achieve convergence with the beam theory result. I am now 
>> using a beam that is 1 x 2 x 20, and using the subdivided rectangle helper 
>> function, I set the number of subdivisions to be 5 and 10 for the width and 
>> height, respectively. I then varied the number of subdivision in the length 
>> between 10 and 1000. The beam theory result is that the shear force has a 
>> magnitude of 0.001 for a displacement on the right side of 0.1. Even at 
>> 1000 subdivisions, the FEM result is 0.00113 (from 0.00129 at 500). The 
>> system has 200 000 degrees of freedom, and the result is still off by 13%. 
>> Is it expected that even to calculate a shear force in this simple problem 
>> that such a large number of degrees of freedom are required?
>>
>>  
>>
>>  
>>
>> Best,
>>
>> Alex
>>
>>  
>>
>> On Monday, June 7, 2021 at 3:39:44 p.m. UTC+2 lian...@gmail.com wrote:
>>
>>
>> Hi Alex,
>>
>> I'm learning deal.ii and trying do the similar verification. If it is 
>> possible for you to share the code with me?
>>
>>  
>>
>> Thank you!
>>
>> Michael
>>
>> On Tuesday, May 11, 2021 at 4:46:55 AM UTC-6 alexanderc...@gmail.com 
>> wrote:
>>
>> Hello,
>>
>>  
>>
>> As a test to validate my code, I am solving the equations for 
>> geometrically nonlinear elasticity (the Saint Venant-Kirchhoff model) for a 
>> beam with a small displacement boundary condition on the right end such 
>> that I can compare with Euler-Bernoulli beam theory. I can compare both the 
>> displacement and the shear force between the FEM solution and the beam 
>> theory solution. In my FEM integration, I output the normal and shear 
>> forces for both sides of the bea

Re: [deal.II] Re: Integrated material and spatial traction forces on boundary not equal

[Wolfgang Bangerth]

On 6/16/21 4:07 PM, Michael Lee wrote:


Is there a way to let deal.II know Trillinos is installed? Do I have to 
reinstall deal.ii?

Or, can I compile the code without using Trillinos with some modification?


No, you have to re-install deal.II so that at compile time, it knows that 
Trilinos exists.


Best
 W.

--

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

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/06a30833-0914-ed13-b4e1-dfd9c2c3b72d%40colostate.edu.

Re: [deal.II] Local Raviart-Thomas space

[Wolfgang Bangerth]

On 6/16/21 2:00 PM, Giselle Sosa Jones wrote:


I am trying to implement a projection of a DG solution to a Raviart-Thomas 
space. For this, I need a local space Q_{p-1,p}(K) \times Q_{p,p-1}(K). I have 
been using FE_DGRaviartThomas of degree p-1, and then I 
use shape_value_component(i, q, 1) when I want to use the x component of the 
test function, and shape_value_component(i,q,0) when I want to use the y 
component of the test function.


I suspect that you made a mistake when you say
  component=1 => x-component
  component=0 => y-component
It should be the other way around.


I have a bug in my code and I am wondering if 
it is coming from the way I am creating this local space. Does this 
construction of the space make sense? Is there maybe a better way of doing it?


It's hard to tell without actually seeing what you do. But I would suggest you 
take a look at step-61, which does the kind of thing you're looking for, I think.


Best
 W.


--

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

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/9fbf3326-5b84-6ca4-15e0-8c33690a9bd5%40colostate.edu.

Re: [deal.II] Cahn-Hilliard Problem

[Wolfgang Bangerth]

On 6/16/21 12:21 PM, Kubra Karayagiz (Alumni) wrote:


Could you please let me know if my understanding of the usage of this function 
is correct? (Note: the use of compute_nl_term() is necessary for the 
implicit-explicit implementation, due to the nonlinear functions)


Kubra,
much debugging will be necessary to determine where the issue lies. I don't 
think any of us will be able to tell just by looking at the code. Have you 
tried a simple case first, say where you start with c=constant and see what 
happens there?


Separately, the problem you are trying to solve has fourth spatial 
derivatives. How do you deal with this?


Best
 W.


--

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

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/3391ca72-d1ca-dbb5-f35a-05185bd488aa%40colostate.edu.

Re: [deal.II] Neumann data and Kelly estimator: can't pass std::map argument

[Wolfgang Bangerth]

On 6/16/21 4:29 AM, bob bill wrote:


*The boundary indicator for Neumann data is 1, and I gave a pointer to a 
constant function as argument. **Unfortunately, the compiler says: "*no 
matching constructor for initialization of 'std::mapFunction<2> *>' (aka 'map *>') 
std::map *>(1,&neumann_func),*"*


*What am I missing?*


It has nothing to do with the function you're trying to call. You just can't 
create a map by saying

  std::map(u,v)
where u is a key and v is a value. But you can
  std::map({{u,v}})

Here, the outer braces enclose a list of pairs, which here has only one 
element and that element is initialized by the pair u,v.


Best
 W.

--

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

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/65597fa1-82de-da60-5379-88137fd14761%40colostate.edu.

Re: [deal.II] Cahn-Hilliard Problem

[Kubra Karayagiz (Alumni)]

 Prof. Bangerth,

Thank you for the prompt response.

On Wed, Jun 16, 2021 at 7:33 PM Wolfgang Bangerth 
wrote:

> On 6/16/21 12:21 PM, Kubra Karayagiz (Alumni) wrote:
> >
> > Could you please let me know if my understanding of the usage of this
> function
> > is correct? (Note: the use of compute_nl_term() is necessary for the
> > implicit-explicit implementation, due to the nonlinear functions)
>
> Kubra,
> much debugging will be necessary to determine where the issue lies. I
> don't
> think any of us will be able to tell just by looking at the code.


I totally understand that. I'll try to figure it out by myself.

Have you
> tried a simple case first, say where you start with c=constant and see
> what
> happens there?
>

I am sorry, I can't quite understand how it could be possible a case with
c=constant, as it is the concentration variable that I am solving for. But,
as I mentioned already my code works properly if I compute the function
fcV(c) within a for loop.  (Based on comparisons with a benchmark study
written in Prisms-PF software adopting deal.ii libraries).

>
> Separately, the problem you are trying to solve has fourth spatial
> derivatives. How do you deal with this?
>

To prevent fourth-order, I split the Cahn-Hilliard equation into two
(Equations (1) and (2) below). I first solve for mu (chemical potential),
then use its solution while solving for c (concentration).

[image: image.png]

Thanks,
Kubra


> Best
>   W.
>
>
> --
> 
> Wolfgang Bangerth  email: bange...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>
> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to the Google Groups
> "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to dealii+unsubscr...@googlegroups.com.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/3391ca72-d1ca-dbb5-f35a-05185bd488aa%40colostate.edu
> .
>

-- 
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/CANnYN%2B06RbCAp3AfeK4jVZweWjH1pdWsondAg9zKuOOiueszUw%40mail.gmail.com.

[deal.II] Extracting a component of FESystem for use elsewhere

[Corbin Foucart]

Hi everyone,

This might be a simple question, but I haven't found a straightforward 
answer in the "handling vector-valued problems" documentation.

If I have a Vector which stores the solution to an FESystem, say 3 
FE_DGQ fields for u, v, p, and I want to extract a new Vector which 
corresponds to only the pressure, what's the best way to do this?

I'm not looking to write the data out for visualization, how to use the 
DataComponentInterpretation classes is clear to me. What I want to do is 
extract a separate Vector which I can then add to another DG scalar 
field of the same type as p in the (u,v,p) system.

I don't need p at the quadrature points, only the nodal values to add to 
another nodal field.

Any advice would be appreciated!
Corbin

-- 
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/6d96604a-cd77-4532-aa16-a6fb638a6ff9n%40googlegroups.com.

[deal.II] Re: Extracting a component of FESystem for use elsewhere

[Corbin Foucart]

I've found one way to do this, but I imagine it's not optimal:

   - call DoFRenumbering::component_wise(fesys_dofh)
   - Do the finite element solve
   - then copy from one vector to another using the number of dofs in that 
   component and an offset
   - This isn't great, however, because I'd rather store the fesys dofh 
   cell-by-cell for memory-friendliness, as almost all other operations are 
   elemental

If there's a better way, I'd be glad to know :)

Best,
Corbin 

On Wednesday, June 16, 2021 at 11:21:56 PM UTC-4 Corbin Foucart wrote:

> Hi everyone,
>
> This might be a simple question, but I haven't found a straightforward 
> answer in the "handling vector-valued problems" documentation.
>
> If I have a Vector which stores the solution to an FESystem, say 3 
> FE_DGQ fields for u, v, p, and I want to extract a new Vector which 
> corresponds to only the pressure, what's the best way to do this?
>
> I'm not looking to write the data out for visualization, how to use the 
> DataComponentInterpretation classes is clear to me. What I want to do is 
> extract a separate Vector which I can then add to another DG scalar 
> field of the same type as p in the (u,v,p) system.
>
> I don't need p at the quadrature points, only the nodal values to add to 
> another nodal field.
>
> Any advice would be appreciated!
> Corbin
>

-- 
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/2bb1442a-830e-46a2-8639-8eec8b101e15n%40googlegroups.com.

Re: [deal.II] Re: Extracting a component of FESystem for use elsewhere

[Paras Kumar]

Hi Corbin,

The cell-wise solution corresponding to the p field could be determined as
follows:
```
dealii::FEValuesExtractorScalar pressureExtractor(2); //assuming you have
three scalar fields and pressure is the last one
dealii::Vector solVecPerCell(nDOFsPerCell);
feValues[pressureExtractor].get_function_values(totalSolution,
solVecPerCell); // here only the values corresponding to pressure are
filled.
```
You could then use this cell-wise values of the p-field in whatever manner
you wish.

Another alternative approach, which I employ, is to use BlockVector object,
wherein one can easily get the total solution vector corresponding to any
of the fields.

Best regards,
Paras

On Thu, Jun 17, 2021 at 7:34 AM Corbin Foucart 
wrote:

> I've found one way to do this, but I imagine it's not optimal:
>
>- call DoFRenumbering::component_wise(fesys_dofh)
>- Do the finite element solve
>- then copy from one vector to another using the number of dofs in
>that component and an offset
>- This isn't great, however, because I'd rather store the fesys dofh
>cell-by-cell for memory-friendliness, as almost all other operations are
>elemental
>
> If there's a better way, I'd be glad to know :)
>
> Best,
> Corbin
>
> On Wednesday, June 16, 2021 at 11:21:56 PM UTC-4 Corbin Foucart wrote:
>
>> Hi everyone,
>>
>> This might be a simple question, but I haven't found a straightforward
>> answer in the "handling vector-valued problems" documentation.
>>
>> If I have a Vector which stores the solution to an FESystem, say
>> 3 FE_DGQ fields for u, v, p, and I want to extract a new Vector
>> which corresponds to only the pressure, what's the best way to do this?
>>
>> I'm not looking to write the data out for visualization, how to use the
>> DataComponentInterpretation classes is clear to me. What I want to do is
>> extract a separate Vector which I can then add to another DG scalar
>> field of the same type as p in the (u,v,p) system.
>>
>> I don't need p at the quadrature points, only the nodal values to add to
>> another nodal field.
>>
>> Any advice would be appreciated!
>> Corbin
>>
> --
> 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/2bb1442a-830e-46a2-8639-8eec8b101e15n%40googlegroups.com
> 
> .
>

-- 
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/CAEU6zmQuPz4OQG4g%2BW%2B20Ax9jEvZzi0fDPb8yKJM391po7h%2BtA%40mail.gmail.com.