From 5c31fcfb5ca6a5cd7e0778a11a34d85858b206c8 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Thu, 5 Aug 2021 18:55:05 +0200
Subject: [PATCH] Add missing dot product for discrete P0 vector functions
---
...EmbeddedIDiscreteFunctionMathFunctions.cpp | 79 ++++++++++++-------
src/scheme/DiscreteFunctionP0Vector.hpp | 14 ++++
2 files changed, 64 insertions(+), 29 deletions(-)
diff --git a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
index a688acf13..d4d10c6c7 100644
--- a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
+++ b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
@@ -3,6 +3,7 @@
#include <language/utils/EmbeddedIDiscreteFunctionUtils.hpp>
#include <mesh/IMesh.hpp>
#include <scheme/DiscreteFunctionP0.hpp>
+#include <scheme/DiscreteFunctionP0Vector.hpp>
#include <scheme/DiscreteFunctionUtils.hpp>
#include <scheme/IDiscreteFunction.hpp>
#include <scheme/IDiscreteFunctionDescriptor.hpp>
@@ -318,43 +319,63 @@ template <size_t Dimension>
std::shared_ptr<const IDiscreteFunction>
dot(const std::shared_ptr<const IDiscreteFunction>& f, const std::shared_ptr<const IDiscreteFunction>& g)
{
- Assert((f->dataType() == ASTNodeDataType::vector_t and f->descriptor().type() == DiscreteFunctionType::P0) and
- (g->dataType() == ASTNodeDataType::vector_t and g->descriptor().type() == DiscreteFunctionType::P0) and
- (f->dataType().dimension() == g->dataType().dimension()));
+ Assert(((f->descriptor().type() == DiscreteFunctionType::P0Vector) and
+ (g->descriptor().type() == DiscreteFunctionType::P0Vector)) or
+ ((f->dataType() == ASTNodeDataType::vector_t and f->descriptor().type() == DiscreteFunctionType::P0) and
+ (g->dataType() == ASTNodeDataType::vector_t and g->descriptor().type() == DiscreteFunctionType::P0) and
+ (f->dataType().dimension() == g->dataType().dimension())));
- using DiscreteFunctionResultType = DiscreteFunctionP0<Dimension, double>;
+ if ((f->descriptor().type() == DiscreteFunctionType::P0Vector) and
+ (g->descriptor().type() == DiscreteFunctionType::P0Vector)) {
+ using DiscreteFunctionResultType = DiscreteFunctionP0<Dimension, double>;
+ using DiscreteFunctionType = DiscreteFunctionP0Vector<Dimension, double>;
- switch (f->dataType().dimension()) {
- case 1: {
- using Rd = TinyVector<1>;
- using DiscreteFunctionType = DiscreteFunctionP0<Dimension, Rd>;
- return std::make_shared<const DiscreteFunctionResultType>(
- dot(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
- }
- case 2: {
- using Rd = TinyVector<2>;
- using DiscreteFunctionType = DiscreteFunctionP0<Dimension, Rd>;
- return std::make_shared<const DiscreteFunctionResultType>(
- dot(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
- }
- case 3: {
- using Rd = TinyVector<3>;
- using DiscreteFunctionType = DiscreteFunctionP0<Dimension, Rd>;
- return std::make_shared<const DiscreteFunctionResultType>(
- dot(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
- }
- default: {
- throw UnexpectedError("invalid data dimension " + operand_type_name(f));
- }
+ const DiscreteFunctionType& f_vector = dynamic_cast<const DiscreteFunctionType&>(*f);
+ const DiscreteFunctionType& g_vector = dynamic_cast<const DiscreteFunctionType&>(*g);
+
+ if (f_vector.size() != g_vector.size()) {
+ throw NormalError("operands have different dimension");
+ } else {
+ return std::make_shared<const DiscreteFunctionResultType>(dot(f_vector, g_vector));
+ }
+
+ } else {
+ using DiscreteFunctionResultType = DiscreteFunctionP0<Dimension, double>;
+
+ switch (f->dataType().dimension()) {
+ case 1: {
+ using Rd = TinyVector<1>;
+ using DiscreteFunctionType = DiscreteFunctionP0<Dimension, Rd>;
+ return std::make_shared<const DiscreteFunctionResultType>(
+ dot(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
+ }
+ case 2: {
+ using Rd = TinyVector<2>;
+ using DiscreteFunctionType = DiscreteFunctionP0<Dimension, Rd>;
+ return std::make_shared<const DiscreteFunctionResultType>(
+ dot(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
+ }
+ case 3: {
+ using Rd = TinyVector<3>;
+ using DiscreteFunctionType = DiscreteFunctionP0<Dimension, Rd>;
+ return std::make_shared<const DiscreteFunctionResultType>(
+ dot(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
+ }
+ default: {
+ throw UnexpectedError("invalid data dimension " + operand_type_name(f));
+ }
+ }
}
}
std::shared_ptr<const IDiscreteFunction>
dot(const std::shared_ptr<const IDiscreteFunction>& f, const std::shared_ptr<const IDiscreteFunction>& g)
{
- if ((f->dataType() == ASTNodeDataType::vector_t and f->descriptor().type() == DiscreteFunctionType::P0) and
- (g->dataType() == ASTNodeDataType::vector_t and g->descriptor().type() == DiscreteFunctionType::P0) and
- (f->dataType().dimension() == g->dataType().dimension())) {
+ if (((f->descriptor().type() == DiscreteFunctionType::P0Vector) and
+ (g->descriptor().type() == DiscreteFunctionType::P0Vector)) or
+ ((f->dataType() == ASTNodeDataType::vector_t and f->descriptor().type() == DiscreteFunctionType::P0) and
+ (g->dataType() == ASTNodeDataType::vector_t and g->descriptor().type() == DiscreteFunctionType::P0) and
+ (f->dataType().dimension() == g->dataType().dimension()))) {
std::shared_ptr mesh = getCommonMesh({f, g});
if (mesh.use_count() == 0) {
diff --git a/src/scheme/DiscreteFunctionP0Vector.hpp b/src/scheme/DiscreteFunctionP0Vector.hpp
index 103718b12..5f00d5964 100644
--- a/src/scheme/DiscreteFunctionP0Vector.hpp
+++ b/src/scheme/DiscreteFunctionP0Vector.hpp
@@ -3,6 +3,7 @@
#include <scheme/IDiscreteFunction.hpp>
+#include <algebra/Vector.hpp>
#include <mesh/Connectivity.hpp>
#include <mesh/ItemArray.hpp>
#include <mesh/Mesh.hpp>
@@ -194,6 +195,19 @@ class DiscreteFunctionP0Vector : public IDiscreteFunction
return product;
}
+ PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
+ dot(const DiscreteFunctionP0Vector& f, const DiscreteFunctionP0Vector& g)
+ {
+ Assert(f.mesh() == g.mesh(), "functions are nor defined on the same mesh");
+ Assert(f.size() == g.size());
+ DiscreteFunctionP0<Dimension, double> result{f.m_mesh};
+ parallel_for(
+ f.m_mesh->numberOfCells(),
+ PUGS_LAMBDA(CellId cell_id) { result[cell_id] = dot(Vector{f[cell_id]}, Vector{g[cell_id]}); });
+
+ return result;
+ }
+
DiscreteFunctionP0Vector(const std::shared_ptr<const MeshType>& mesh, size_t size)
: m_mesh{mesh}, m_cell_arrays{mesh->connectivity(), size}
{}
--
GitLab