-
Notifications
You must be signed in to change notification settings - Fork 3
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Fe br #8
base: master
Are you sure you want to change the base?
Fe br #8
Conversation
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Looks overall pretty good. Just a few comments. Would you like to extract the parts that extend the mapping_type to a vector and submit that individually? These things are often easier to get through patch review in smaller increments.
include/deal.II/fe/fe_poly_tensor.h
Outdated
MappingType mapping_type = this->mapping_type[i]; | ||
|
||
if (mapping_type != mapping_none) | ||
update_transformed_shape_values = true; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
In essence, what you're doing here is determine whether any of the elements of the vector are different from mapping_none
. (And similarly in the next few lines.) You can write this more concisely in the following way (see https://en.cppreference.com/w/cpp/algorithm/find):
const book update_transformed_shape_values = (std::find_if_not (this->mapping_type.begin(),
this->mapping_type.begin(),
[](const MappingType t) { return t != mapping_none; })
!= this->mapping_type.end());
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This command is a mouthful. Let me see... so this command locates any entry between begin or end satisfies the inline function (what's the C++ word for that?) that returns true if the mapping type is not mapping_none. I assume the default return is .end()
if nothing is found?
const bool update_transformed_shape_values = (std::find_if (this->mapping_type.begin(),
this->mapping_type.end(),
[](const MappingType t) { return t != mapping_none; })
!= this->mapping_type.end());
Should that be .end()
on the second line? Also, if we use return t != mapping_none
, shouldn't we use find_if instead of find_if_not?
And if I were to do this for transformed_shape_grads, would it look like this?
const bool update_transformed_shape_values = (std::find_if (this->mapping_type.begin(),
this->mapping_type.end(),
[](const MappingType t) { return (t == mapping_raviart_thomas ||
t == mapping_piola ||
t == mapping_nedelec ||
t == mapping_contravariant); })
!= this->mapping_type.end());
If so, then I think I understand how this works and I can implement it.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Oh is that called a "lambda expression"? This is something I haven't explored much of before.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes:
- That's a lambda function -- a function without a name defined in the place where you need it.
- And yes, correct, it should have been
.end()
in the second line. - And yes, if
find_if
can't find an entry that satisfies the lambda function, then it returns.end()
.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
include/deal.II/fe/fe_poly_tensor.h
Outdated
(mapping_type == mapping_piola) || | ||
(mapping_type == mapping_nedelec) || | ||
(mapping_type == mapping_contravariant)) | ||
update_transformed_shape_grads = true; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This would work with find_if (...)
.
include/deal.II/fe/fe_poly_tensor.h
Outdated
update_transformed_shape_grads = true; | ||
|
||
if (mapping_type != mapping_none) | ||
update_transformed_shape_hessian_tensors = true; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is the exact same condition as you had above for update_transformed_shape_values
, so you will get the exact same answer. Is this on purpose?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Ah it was not on purpose. Seeing the original command looked like
if (mapping_type != mapping_none)
data->untransformed_shape_hessian_tensors.resize(n_q_points);
I can just absorb that boolean as being equal to the update_transformed_shape_values
boolean since both only check if the mapping is not mapping_none
.
source/fe/fe_bernardi_raugel.cc
Outdated
mapping_none, mapping_none, | ||
mapping_none, mapping_none, | ||
mapping_piola, mapping_piola, | ||
mapping_piola, mapping_piola}; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I believe you'll have to update that, right?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes, I will revert it to {mapping_none}
once the fill_fe_values overload is finished.
source/fe/fe_bernardi_raugel.cc
Outdated
mapping_none, mapping_none, mapping_none, | ||
mapping_piola, mapping_piola, | ||
mapping_piola, mapping_piola, | ||
mapping_piola, mapping_piola}; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
add the following:
else
Assert (false, ExcNotImplemented());
to make sure that if the code is ever used for anything other than 2d and 3d, you will be notified that that case has not been implemented.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
On line 53 I have Assert(dim == 2 || dim == 3, ExcImpossibleInDim(dim));
Should I still add this one as well?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes, I think that would be useful. The assertion at the top is an expression of the general idea that this function isn't implemented for certain cases, whereas adding the assertion here would point to a specific place where could would need to be augmented if one wanted to extend the function.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Okay got it!
} | ||
|
||
|
||
|
||
// template<int spacedim> |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think you could do it with this template argument as well.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I only hesitated because on line 73 there's a const unsigned int spacedim = 2;
Should I still throw that template in there?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Just remove that line. I don't think that the function should behave differently based on the space dimension; only the dimension should matter.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Okay I updated that. Would I need to update any other files since I'm turning this into a templated function? I know there are .inst files and they relate with templates but I'm unsure how those things work.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
No, the function is used only in this one file. The compiler can figure this out. It only gets complicated when you declare a template in one (.h) file, define (i.e., implement) it in a .cc file, and want to use it in another .cc file.
} | ||
} | ||
} | ||
} | ||
|
||
|
||
|
||
template <int spacedim> |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
yes, just like you do here
source/fe/fe_poly_tensor.cc
Outdated
face_sign[f * dofs_per_face + j] = -1.0; | ||
if (mapping_type.size() > 1 ? | ||
mapping_type[cell_j] == mapping_raviart_thomas : | ||
mapping_type[0] == mapping_raviart_thomas) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It's not wrong to use the ?:
operator in an if
clause, but it's often easier to read to write things out:
if ((cond && expr1) || (!cond && expr2))
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Or here maybe:
if ((mapping_type.size() > 1 ? mapping_type[cell_j] : mapping_type[0]) == mapping_raviart_thomas)
This probably expresses most clearly what you are doing.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes I like that second option.
this->mapping_type, | ||
fe_data.sign_change); | ||
|
||
internal::FE_PolyTensor::get_face_sign_change_nedelec(cell, |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is clearer now. I wonder whether the two functions could be merged into one, though I'd leave this for a follow-up patch.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I thought of that when I saw these side-by-side. If these were merged into one function, it would likely need to be a function that handles rt, nedelec, and is easy to modify in case future elements require something similar.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes. Leave it as is for the moment. We can always come back to it.
source/fe/fe_poly_tensor.cc
Outdated
|
||
for (unsigned int i = 0; i < this->dofs_per_cell; ++i) | ||
{ | ||
MappingType mapping_type = |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
here and below, mark this variable as const
since it won't change in this block.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Got it!
Apologies for taking this long to get to looking at your changes. Let me know if you want me to read through them again! |
Fixes: ``` 134: ==1438450==ERROR: AddressSanitizer: stack-use-after-return on address 0x7fa4f4d4bfb0 at pc 0x56248d8610c2 bp 0x7ffeebad19b0 sp 0x7ffeebad19a8 134: READ of size 1 at 0x7fa4f4d4bfb0 thread T0 134: #0 0x56248d8610c1 in ExponentialDecay::ExponentialDecay(double, dealii::PETScWrappers::TimeStepperData const&, bool, bool, int)::'lambda'(double, dealii::[148/1839] ers::MPI::Vector const&, dealii::PETScWrappers::MPI::Vector const&, double, dealii::PETScWrappers::MatrixBase&, dealii::PETScWrappers::MatrixBase&)::operator()(double, deali i::PETScWrappers::MPI::Vector const&, dealii::PETScWrappers::MPI::Vector const&, double, dealii::PETScWrappers::MatrixBase&, dealii::PETScWrappers::MatrixBase&) const /srv/t estsuite/dealii/tests/petsc/petsc_ts_05.cc:91:15 134: #1 0x7fa538cf812b in std::__1::__function::__value_func<void (double, dealii::PETScWrappers::MPI::Vector const&, dealii::PETScWrappers::MPI::Vector const&, double, dealii::PETScWrappers::MatrixBase&, dealii::PETScWrappers::MatrixBase&)>::operator()[abi:v160006](double&&, dealii::PETScWrappers::MPI::Vector const&, dealii::PETScWrappers: :MPI::Vector const&, double&&, dealii::PETScWrappers::MatrixBase&, dealii::PETScWrappers::MatrixBase&) const /usr/include/c++/v1/__functional/function.h:510:16 134: #2 0x7fa538cf812b in std::__1::function<void (double, dealii::PETScWrappers::MPI::Vector const&, dealii::PETScWrappers::MPI::Vector const&, double, dealii::PETScWra ppers::MatrixBase&, dealii::PETScWrappers::MatrixBase&)>::operator()(double, dealii::PETScWrappers::MPI::Vector const&, dealii::PETScWrappers::MPI::Vector const&, double, de alii::PETScWrappers::MatrixBase&, dealii::PETScWrappers::MatrixBase&) const /usr/include/c++/v1/__functional/function.h:1156:12 134: #3 0x7fa538cf812b in int dealii::PETScWrappers::call_and_possibly_capture_ts_exception<std::__1::function<void (double, dealii::PETScWrappers::MPI::Vector const&, d ealii::PETScWrappers::MPI::Vector const&, double, dealii::PETScWrappers::MatrixBase&, dealii::PETScWrappers::MatrixBase&)>, double&, dealii::PETScWrappers::MPI::Vector&, dea lii::PETScWrappers::MPI::Vector&, double&, dealii::PETScWrappers::MatrixBase&, dealii::PETScWrappers::MatrixBase&>(std::__1::function<void (double, dealii::PETScWrappers::MP I::Vector const&, dealii::PETScWrappers::MPI::Vector const&, double, dealii::PETScWrappers::MatrixBase&, dealii::PETScWrappers::MatrixBase&)> const&, std::exception_ptr&, st d::__1::function<void ()> const&, double&, dealii::PETScWrappers::MPI::Vector&, dealii::PETScWrappers::MPI::Vector&, double&, dealii::PETScWrappers::MatrixBase&, dealii::PET ScWrappers::MatrixBase&) /srv/testsuite/dealii/include/deal.II/lac/petsc_ts.templates.h:87:9 134: #4 0x7fa538cf78cc in dealii::PETScWrappers::TimeStepper<dealii::PETScWrappers::MPI::Vector, dealii::PETScWrappers::MatrixBase, dealii::PETScWrappers::MatrixBase>::s olve(dealii::PETScWrappers::MPI::Vector&)::'lambda'(_p_TS*, double, _p_Vec*, _p_Vec*, double, _p_Mat*, _p_Mat*, void*)::operator()(_p_TS*, double, _p_Vec*, _p_Vec*, double, _p_Mat*, _p_Mat*, void*) const /srv/testsuite/dealii/include/deal.II/lac/petsc_ts.templates.h:489:26 134: #5 0x7fa538cf73d7 in dealii::PETScWrappers::TimeStepper<dealii::PETScWrappers::MPI::Vector, dealii::PETScWrappers::MatrixBase, dealii::PETScWrappers::MatrixBase>::s olve(dealii::PETScWrappers::MPI::Vector&)::'lambda'(_p_TS*, double, _p_Vec*, _p_Vec*, double, _p_Mat*, _p_Mat*, void*)::__invoke(_p_TS*, double, _p_Vec*, _p_Vec*, double, _p _Mat*, _p_Mat*, void*) /srv/testsuite/dealii/include/deal.II/lac/petsc_ts.templates.h:472:31 134: #6 0x7fa504afbcac in TSComputeIJacobian (/usr/lib64/libpetsc.so.3.19+0x10fbcac) 134: #7 0x7fa504aa1672 (/usr/lib64/libpetsc.so.3.19+0x10a1672) 134: #8 0x7fa504afd307 in SNESTSFormJacobian (/usr/lib64/libpetsc.so.3.19+0x10fd307) 134: #9 0x7fa504a43e66 in SNESComputeJacobian (/usr/lib64/libpetsc.so.3.19+0x1043e66) 134: #10 0x7fa504a35e98 in SNESSolve_NEWTONLS (/usr/lib64/libpetsc.so.3.19+0x1035e98) 134: #11 0x7fa504a4968e in SNESSolve (/usr/lib64/libpetsc.so.3.19+0x104968e) 134: #12 0x7fa504aa25e4 (/usr/lib64/libpetsc.so.3.19+0x10a25e4) 134: #13 0x7fa504aa014a (/usr/lib64/libpetsc.so.3.19+0x10a014a) 134: #14 0x7fa504b03013 in TSStep (/usr/lib64/libpetsc.so.3.19+0x1103013) 134: #15 0x7fa504b03e02 in TSSolve (/usr/lib64/libpetsc.so.3.19+0x1103e02) 134: #16 0x7fa538cd2ccc in dealii::PETScWrappers::TimeStepper<dealii::PETScWrappers::MPI::Vector, dealii::PETScWrappers::MatrixBase, dealii::PETScWrappers::MatrixBase>:: solve(dealii::PETScWrappers::MPI::Vector&) /srv/testsuite/dealii/include/deal.II/lac/petsc_ts.templates.h:880:24 134: #17 0x56248d85849f in ExponentialDecay::run() /srv/testsuite/dealii/tests/petsc/petsc_ts_05.cc:182:32 134: #18 0x56248d855ca6 in main /srv/testsuite/dealii/tests/petsc/petsc_ts_05.cc:242:20 134: #19 0x7fa502a50989 (/usr/lib64/libc.so.6+0x23989) 134: #20 0x7fa502a50a44 in __libc_start_main (/usr/lib64/libc.so.6+0x23a44) 134: #21 0x56248d73e2d0 in _start (/srv/temp/build/tests/petsc/petsc_ts_05.debug/petsc_ts_05.debug+0x742d0) ```
``` ==1080297==ERROR: AddressSanitizer: attempting free on address which was not malloc()-ed: 0x60200003da70 in thread T0 #0 0x55bcdd907b7d in operator delete(void*) (/srv/temp/testsuite-IQZ1b8kK/build/tests/scalapack/scalapack_06b.debug/scalapack_06b.debug+0x17cb7d) #1 0x7fc52a4d2047 in void std::__1::__libcpp_operator_delete[abi:v160006]<void*>(void*) /usr/include/c++/v1/new:276:3 #2 0x7fc52a4d2047 in void std::__1::__do_deallocate_handle_size[abi:v160006]<>(void*, unsigned long) /usr/include/c++/v1/new:300:10 #3 0x7fc52a4d2047 in std::__1::__libcpp_deallocate[abi:v160006](void*, unsigned long, unsigned long) /usr/include/c++/v1/new:316:14 #4 0x7fc52a4d2047 in std::__1::allocator<int>::deallocate[abi:v160006](int*, unsigned long) /usr/include/c++/v1/__memory/allocator.h:131:13 #5 0x7fc52a4d2047 in std::__1::allocator_traits<std::__1::allocator<int>>::deallocate[abi:v160006](std::__1::allocator<int>&, int*, unsigned long) /usr/include/c++/v1/__memory/allocator_traits.h:288:13 #6 0x7fc52a4d2047 in std::__1::__split_buffer<int, std::__1::allocator<int>&>::~__split_buffer() /usr/include/c++/v1/__split_buffer:362:9 #7 0x7fc52a4d2047 in std::__1::vector<int, std::__1::allocator<int>>::__append(unsigned long) /usr/include/c++/v1/vector:1049:5 #8 0x7fc52a8282d2 in std::__1::vector<int, std::__1::allocator<int>>::resize(unsigned long) /usr/include/c++/v1/vector:1910:15 #9 0x7fc52a8282d2 in dealii::ScaLAPACKMatrix<double>::eigenpairs_symmetric(bool, std::__1::pair<unsigned int, unsigned int> const&, std::__1::pair<double, double> const&) /srv/temp/testsuite-IQZ1b8kK/dealii/source/lac/scalapack.cc:1684:17 #10 0x7fc52a826813 in dealii::ScaLAPACKMatrix<double>::eigenpairs_symmetric_by_index(std::__1::pair<unsigned int, unsigned int> const&, bool) /srv/temp/testsuite-IQZ1b8kK/dealii/source/lac/scalapack.cc:1446:12 #11 0x55bcdd911907 in void test<double>(unsigned int, unsigned int, double) /srv/temp/testsuite-IQZ1b8kK/dealii/tests/scalapack/scalapack_06b.cc:134:21 #12 0x55bcdd90ebc8 in main /srv/temp/testsuite-IQZ1b8kK/dealii/tests/scalapack/scalapack_06b.cc:207:9 #13 0x7fc4f5250989 (/usr/lib64/libc.so.6+0x23989) #14 0x7fc4f5250a44 in __libc_start_main (/usr/lib64/libc.so.6+0x23a44) #15 0x55bcdd7f73f0 in _start (/srv/temp/testsuite-IQZ1b8kK/build/tests/scalapack/scalapack_06b.debug/scalapack_06b.debug+0x6c3f0) ```
``` ==1080297==ERROR: AddressSanitizer: attempting free on address which was not malloc()-ed: 0x60200003da70 in thread T0 #0 0x55bcdd907b7d in operator delete(void*) (/srv/temp/testsuite-IQZ1b8kK/build/tests/scalapack/scalapack_06b.debug/scalapack_06b.debug+0x17cb7d) #1 0x7fc52a4d2047 in void std::__1::__libcpp_operator_delete[abi:v160006]<void*>(void*) /usr/include/c++/v1/new:276:3 #2 0x7fc52a4d2047 in void std::__1::__do_deallocate_handle_size[abi:v160006]<>(void*, unsigned long) /usr/include/c++/v1/new:300:10 #3 0x7fc52a4d2047 in std::__1::__libcpp_deallocate[abi:v160006](void*, unsigned long, unsigned long) /usr/include/c++/v1/new:316:14 #4 0x7fc52a4d2047 in std::__1::allocator<int>::deallocate[abi:v160006](int*, unsigned long) /usr/include/c++/v1/__memory/allocator.h:131:13 #5 0x7fc52a4d2047 in std::__1::allocator_traits<std::__1::allocator<int>>::deallocate[abi:v160006](std::__1::allocator<int>&, int*, unsigned long) /usr/include/c++/v1/__memory/allocator_traits.h:288:13 #6 0x7fc52a4d2047 in std::__1::__split_buffer<int, std::__1::allocator<int>&>::~__split_buffer() /usr/include/c++/v1/__split_buffer:362:9 #7 0x7fc52a4d2047 in std::__1::vector<int, std::__1::allocator<int>>::__append(unsigned long) /usr/include/c++/v1/vector:1049:5 #8 0x7fc52a8282d2 in std::__1::vector<int, std::__1::allocator<int>>::resize(unsigned long) /usr/include/c++/v1/vector:1910:15 #9 0x7fc52a8282d2 in dealii::ScaLAPACKMatrix<double>::eigenpairs_symmetric(bool, std::__1::pair<unsigned int, unsigned int> const&, std::__1::pair<double, double> const&) /srv/temp/testsuite-IQZ1b8kK/dealii/source/lac/scalapack.cc:1684:17 #10 0x7fc52a826813 in dealii::ScaLAPACKMatrix<double>::eigenpairs_symmetric_by_index(std::__1::pair<unsigned int, unsigned int> const&, bool) /srv/temp/testsuite-IQZ1b8kK/dealii/source/lac/scalapack.cc:1446:12 #11 0x55bcdd911907 in void test<double>(unsigned int, unsigned int, double) /srv/temp/testsuite-IQZ1b8kK/dealii/tests/scalapack/scalapack_06b.cc:134:21 #12 0x55bcdd90ebc8 in main /srv/temp/testsuite-IQZ1b8kK/dealii/tests/scalapack/scalapack_06b.cc:207:9 #13 0x7fc4f5250989 (/usr/lib64/libc.so.6+0x23989) #14 0x7fc4f5250a44 in __libc_start_main (/usr/lib64/libc.so.6+0x23a44) #15 0x55bcdd7f73f0 in _start (/srv/temp/testsuite-IQZ1b8kK/build/tests/scalapack/scalapack_06b.debug/scalapack_06b.debug+0x6c3f0) ```
No description provided.