I found a fix for this by editing the deal.ii source code. I am very new to
deal.ii so it is extremely likely that there is a way to accomplish what I
want without touching the deal.ii source or there might be a better way of
editing the source. So, here is the fix that worked for me. Let me know if
there is a better solution.
I first implemented a new reinit function for IntegrationInfo:
/**
* Reinitialize internal data structures for use on a cell.
*/
template <typename number, class ITERATOR>
void
reinit(ITERATOR cell, const DoFInfo<dim, spacedim, number> &info);
template <int dim, int spacedim>
template <typename number, class ITERATOR>
inline void
IntegrationInfo<dim, spacedim>::reinit(ITERATOR cell,
const DoFInfo<dim, spacedim, number> &info)
{
for (unsigned int i = 0; i < fevalv.size(); ++i)
{
FEValuesBase<dim, spacedim> &febase = *fevalv[i];
if (info.sub_number != numbers::invalid_unsigned_int)
{
// This is a subface
FESubfaceValues<dim, spacedim> &fe =
dynamic_cast<FESubfaceValues<dim, spacedim> &>(febase);
fe.reinit(cell, info.face_number, info.sub_number);
}
else if (info.face_number != numbers::invalid_unsigned_int)
{
// This is a face
FEFaceValues<dim, spacedim> &fe =
dynamic_cast<FEFaceValues<dim, spacedim> &>(febase);
fe.reinit(cell, info.face_number);
}
else
{
// This is a cell
FEValues<dim, spacedim> &fe =
dynamic_cast<FEValues<dim, spacedim> &>(febase);
fe.reinit(cell);
}
}
const bool split_fevalues = info.block_info != nullptr;
if (!global_data->empty())
fill_local_data(info, split_fevalues);
}
And then I changed the MeshWorker::loop function on lines 231-232 of
deal.ii/meshworker/loop.h to
if (integrate_cell)
info.cell.reinit(cell, dof_info.cell);
Cheers,
Andy
On Thursday, January 30, 2020 at 3:39:28 PM UTC-5, Andrew Davis wrote:
>
> Fair point. I reinstalled in Debug mode and now I'm failing at that
> Assert. Here is the error message (using the same code I attached
> previously):
>
> An error occurred in line <2891> of file <dealii/source/fe/fe_values.cc>
> in function
> dealii::types::global_dof_index dealii::FEValuesBase<dim,
> spacedim>::TriaCellIterator::n_dofs_for_dof_handler() const [with int dim =
> 2; int spacedim = 2; dealii::types::global_dof_index = unsigned int]
> The violated condition was:
> false
> Additional information:
> You have previously called the FEValues::reinit function with a
> cell iterator of type Triangulation<dim,spacedim>::cell_iterator. However,
> when you do this, you cannot call some functions in the FEValues
> class, such as the get_function_values/gradients/hessians/third_derivatives
> functions. If you need these functions, then you need to call
> FEValues::reinit with an iterator type that allows to extract
> degrees of freedom, such as DoFHandler<dim,spacedim>::cell_iterator.
>
> I'm going to dig into this further but does anyone know if this is a bug
> or user error? If it is user error does anyone know how to fix this?
>
> Cheers,
> Andy
>
> On Thursday, January 30, 2020 at 2:38:00 PM UTC-5, Bruno Turcksin wrote:
>>
>> On Thursday, January 30, 2020 at 2:28:55 PM UTC-5, Andrew Davis wrote:
>>>
>>>
>>> from the file fe_values.impl.2.inst.in. The Assert(false) seems strange
>>> to me---why doesn't it crash? (I'm in Release mode so that could answer
>>> that). However, should it be calling a different function?
>>>
>>> You should always debug your code in Debug mode. It's very possible that
>> there is an assert that will catch a problem earlier in the code and so the
>> code after that assert makes no sense because you are not supposed to get
>> there.
>>
>> 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/13e914ec-aa34-4f8e-86fc-84b97b74f9f5%40googlegroups.com.
