From bffc0c7fa2e148a898d6e24eb791fa5961d5d556 Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Wed, 16 Sep 2020 18:09:38 -0500 Subject: [PATCH] Remove boundary fix We compute the boundary normal by averaging the normal over the quadrature points weighted with the shape function. It turns out that this is sufficient to preserve mass and total energy. No change of the c_ij's are needed. --- examples/step-69/step-69.cc | 116 ------------------------------------ 1 file changed, 116 deletions(-) diff --git a/examples/step-69/step-69.cc b/examples/step-69/step-69.cc index c6c428642f..f513ab165f 100644 --- a/examples/step-69/step-69.cc +++ b/examples/step-69/step-69.cc @@ -1295,122 +1295,6 @@ namespace Step69 normal /= (normal.norm() + std::numeric_limits::epsilon()); } } - - // As commented in the introduction (see section titled "Stable boundary - // conditions and conservation properties") we use three different types of - // boundary conditions: essential-like boundary conditions (we prescribe a - // state in the left portion of our domain), outflow boundary conditions - // (also called "do-nothing" boundary conditions) in the right portion of - // the domain, and "reflecting" boundary conditions (also called "free - // slip" boundary conditions). With these boundary conditions we should - // not expect any form of conservation to hold. - // - // However, if we were to use reflecting boundary conditions - // $\mathbf{m} \cdot \boldsymbol{\nu}_i =0$ on the entirety of the - // boundary we should preserve the density and total (mechanical) - // energy. This requires us to modify the $\mathbf{c}_{ij}$ vectors at - // the boundary as follows @cite GuermondEtAl2018 : - // - // @f{align*} - // \mathbf{c}_{ij} \, +\!\!= \int_{\partial \Omega} - // (\boldsymbol{\nu}_j - \boldsymbol{\nu}(\mathbf{s})) \phi_i \phi_j \, - // \mathrm{d}\mathbf{s} \ \text{whenever} \ \mathbf{x}_i \text{ and } - // \mathbf{x}_j \text{ lie in the boundary.} - // @f} - // - // The ideas repeat themselves: we use WorkStream in order to compute - // this correction, most of the following code is about the definition - // of the worker local_assemble_system(). - { - TimerOutput::Scope scope(computing_timer, - "offline_data - fix slip boundary c_ij"); - - const auto local_assemble_system = // - [&](const typename DoFHandler::cell_iterator &cell, - MeshWorker::ScratchData & scratch, - CopyData & copy) { - copy.is_artificial = cell->is_artificial(); - - if (copy.is_artificial) - return; - - for (auto &matrix : copy.cell_cij_matrix) - matrix.reinit(dofs_per_cell, dofs_per_cell); - - copy.local_dof_indices.resize(dofs_per_cell); - cell->get_dof_indices(copy.local_dof_indices); - std::transform(copy.local_dof_indices.begin(), - copy.local_dof_indices.end(), - copy.local_dof_indices.begin(), - [&](types::global_dof_index index) { - return partitioner->global_to_local(index); - }); - - for (auto &matrix : copy.cell_cij_matrix) - matrix = 0.; - - for (const auto f : cell->face_indices()) - { - const auto face = cell->face(f); - const auto id = face->boundary_id(); - - if (!face->at_boundary()) - continue; - - if (id != Boundaries::free_slip) - continue; - - const auto &fe_face_values = scratch.reinit(cell, f); - - const unsigned int n_face_q_points = - fe_face_values.get_quadrature().size(); - - for (unsigned int q = 0; q < n_face_q_points; ++q) - { - const auto JxW = fe_face_values.JxW(q); - const auto normal_q = fe_face_values.normal_vector(q); - - for (unsigned int j = 0; j < dofs_per_cell; ++j) - { - if (!discretization->finite_element.has_support_on_face( - j, f)) - continue; - - const auto &normal_j = std::get<0>( - boundary_normal_map[copy.local_dof_indices[j]]); - - const auto value_JxW = - fe_face_values.shape_value(j, q) * JxW; - - for (unsigned int i = 0; i < dofs_per_cell; ++i) - { - const auto value = fe_face_values.shape_value(i, q); - - /* This is the correction of the boundary c_ij */ - for (unsigned int d = 0; d < dim; ++d) - copy.cell_cij_matrix[d](i, j) += - (normal_j[d] - normal_q[d]) * (value * value_JxW); - } /* i */ - } /* j */ - } /* q */ - } /* f */ - }; - - const auto copy_local_to_global = [&](const CopyData ©) { - if (copy.is_artificial) - return; - - for (int k = 0; k < dim; ++k) - cij_matrix[k].add(copy.local_dof_indices, copy.cell_cij_matrix[k]); - }; - - WorkStream::run(dof_handler.begin_active(), - dof_handler.end(), - local_assemble_system, - copy_local_to_global, - scratch_data, - CopyData()); - } } // At this point we are very much done with anything related to offline data. -- 2.39.5