30 std::vector<FluxVectorType>
const& ip_flux_vector,
31 Eigen::MatrixBase<Derived>& laplacian_matrix)
33 auto const& Ns = shape_matrix_cache.
NsHigherOrder<MeshElementType>();
35 for (std::size_t ip = 0; ip < ip_flux_vector.size(); ++ip)
37 auto const& ip_data = ip_data_vector[ip];
38 auto const w = ip_data.integration_weight;
39 auto const& dNdx = ip_data.dNdx;
40 auto const&
N = Ns[ip];
41 laplacian_matrix.noalias() +=
42 (
N.transpose() * ip_flux_vector[ip].transpose() * dNdx * w)
43#if defined(__GNUC__) && (__GNUC__ == 14) && (__GNUC_MINOR__ == 2) && \
44 (__GNUC_PATCHLEVEL__ == 1)
53 std::vector<FluxVectorType>
const& ip_flux_vector,
54 Eigen::MatrixBase<Derived>& laplacian_matrix)
56 for (std::size_t ip = 0; ip < ip_flux_vector.size(); ++ip)
58 auto const& ip_data = ip_data_vector[ip];
59 auto const w = ip_data.integration_weight;
60 auto const&
N = ip_data.N;
61 auto const& dNdx = ip_data.dNdx;
62 laplacian_matrix.noalias() +=
63 N.transpose() * ip_flux_vector[ip].transpose() * dNdx * w;
69 Eigen::MatrixBase<Derived>& laplacian_matrix)
71 Eigen::VectorXd
const down_mask =
72 (quasi_nodal_flux.array() < 0).cast<
double>();
73 Eigen::VectorXd
const down = quasi_nodal_flux.cwiseProduct(down_mask);
75 double const q_in = -down.sum();
76 if (q_in < std::numeric_limits<double>::epsilon())
81 Eigen::VectorXd
const up_mask =
82 (quasi_nodal_flux.array() >= 0).cast<
double>();
83 Eigen::VectorXd
const up = quasi_nodal_flux.cwiseProduct(up_mask);
85 laplacian_matrix.diagonal().noalias() += up;
86 laplacian_matrix.noalias() += down * up.transpose() / q_in;
91 std::vector<FluxVectorType>
const& ip_flux_vector,
92 Eigen::MatrixBase<Derived>& laplacian_matrix)
94 Eigen::VectorXd quasi_nodal_flux =
95 Eigen::VectorXd::Zero(laplacian_matrix.rows());
97 for (std::size_t ip = 0; ip < ip_flux_vector.size(); ++ip)
99 auto const& ip_data = ip_data_vector[ip];
100 auto const w = ip_data.integration_weight;
101 auto const& dNdx = ip_data.dNdx;
103 quasi_nodal_flux.noalias() -= ip_flux_vector[ip].transpose() * dNdx * w;
115 IPData
const& ip_data_vector,
117 std::vector<FluxVectorType>
const& ip_flux_vector,
118 double const average_velocity,
119 Eigen::MatrixBase<Derived>& laplacian_matrix)
122 [&](
auto&& stabilizer)
124 using Stabilizer = std::decay_t<
decltype(stabilizer)>;
125 if constexpr (std::is_same_v<Stabilizer, FullUpwind>)
127 if (average_velocity > stabilizer.getCutoffVelocity())
136 ip_data_vector, shape_matrix_cache, ip_flux_vector,
144 IPData
const& ip_data_vector,
145 std::vector<FluxVectorType>
const& ip_flux_vector,
146 double const average_velocity,
147 Eigen::MatrixBase<Derived>& laplacian_matrix)
150 [&](
auto&& stabilizer)
152 using Stabilizer = std::decay_t<
decltype(stabilizer)>;
153 if constexpr (std::is_same_v<Stabilizer, FullUpwind>)
155 if (average_velocity > stabilizer.getCutoffVelocity())
164 ip_data_vector, ip_flux_vector, laplacian_matrix);