diff --git a/CMakeLists.txt b/CMakeLists.txt index dd9a408..17b7b51 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -133,6 +133,16 @@ if(OpenMP_CXX_FOUND) message(STATUS "Found OpenMP") endif() +if(OpenMP_CXX_FOUND) + message(STATUS "Found OpenMP: ${OpenMP_CXX_FLAGS}") + + add_library(cfem_openmp_interface INTERFACE) + + target_link_libraries(cfem_openmp_interface INTERFACE OpenMP::OpenMP_CXX) + + target_compile_definitions(cfem_openmp_interface INTERFACE WITH_OPENMP) +endif() + # ------------------------------------------------------------ # Modular Libraries & Applications # ------------------------------------------------------------ @@ -156,4 +166,4 @@ if(BUILD_TESTING) include(GoogleTest) add_subdirectory(tests) -endif() \ No newline at end of file +endif() diff --git a/apps/main.cpp b/apps/main.cpp index f1db4ae..3d461c3 100644 --- a/apps/main.cpp +++ b/apps/main.cpp @@ -53,7 +53,7 @@ static char help[] = "CFEM Example: Field projection on a quadratic mesh.\n\n"; int main(int argc, char **argv) { - //omp_set_num_threads(1); + // omp_set_num_threads(4); Environment::initialize(argc, argv, help); { #define CFEM_USE_DYNAMIC_SPACE_CACHE diff --git a/libs/expressions/CMakeLists.txt b/libs/expressions/CMakeLists.txt index d826cd3..695ae9e 100644 --- a/libs/expressions/CMakeLists.txt +++ b/libs/expressions/CMakeLists.txt @@ -35,4 +35,5 @@ target_link_libraries(cfem_expressions ${MPFR_LIB} ${MPC_LIB} # ${LLVM_SYSTEM_LIBS} -) \ No newline at end of file + cfem_openmp_interface +) diff --git a/libs/expressions/include/CfemExpression.h b/libs/expressions/include/CfemExpression.h index 7313b70..c61895a 100644 --- a/libs/expressions/include/CfemExpression.h +++ b/libs/expressions/include/CfemExpression.h @@ -97,7 +97,7 @@ namespace cfem::sym g->evaluate(cellId, xi, out, ctx); return; } - CFEM_ERROR << "Gradient evaluation failed."; + CFEM_THROW("Gradient evaluation failed."); } virtual void prepare() const {}; @@ -128,7 +128,7 @@ namespace cfem::sym m_gradient = createGradient(); } if ( !m_gradient ){ - CFEM_ERROR << "Could not differentiate expression: " << toString(); + CFEM_THROW("Could not differentiate expression: " << toString()); } return m_gradient; } @@ -138,7 +138,7 @@ namespace cfem::sym m_timeDeriv = createTimeDerivative(); } if ( !m_timeDeriv ){ - CFEM_ERROR << "Could not differentiate expression: " << toString(); + CFEM_THROW("Could not differentiate expression: " << toString()); } return m_timeDeriv; } @@ -180,7 +180,7 @@ namespace cfem::sym { // Otherwise, check if it's already used by someone else to avoid collisions if (active_variables.contains(name)) { - CFEM_ERROR << "Variable name '" << name << "' is already in use."; + CFEM_THROW("Variable name '" << name << "' is already in use."); } active_variables.insert(name); } diff --git a/libs/expressions/include/DivergenceExpression.h b/libs/expressions/include/DivergenceExpression.h index 31c8b19..2d0c2b7 100644 --- a/libs/expressions/include/DivergenceExpression.h +++ b/libs/expressions/include/DivergenceExpression.h @@ -50,7 +50,7 @@ namespace cfem::sym } std::shared_ptr createGradient() const override { - CFEM_ERROR << "The gradient of a divergence node requires a dedicated GradExpression fallback."; + CFEM_THROW("The gradient of a divergence node requires a dedicated GradExpression fallback."); return nullptr; } @@ -61,13 +61,13 @@ namespace cfem::sym m_dimIn = m_expr->getDimension(); if (m_dimIn == 1) { - CFEM_ERROR << "Divergence cannot be applied to a scalar field: " << m_expr->toString(); + CFEM_THROW("Divergence cannot be applied to a scalar field: " << m_expr->toString()); } else if (m_dimIn == 3) { m_dimOut = 1; // Vector -> Scalar } else if (m_dimIn == 9) { m_dimOut = 3; // Tensor -> Vector } else { - CFEM_ERROR << "Unsupported dimension for divergence: " << m_dimIn; + CFEM_THROW("Unsupported dimension for divergence: " << m_dimIn); } m_dx = m_expr->diff(0); diff --git a/libs/expressions/include/GradExpression.h b/libs/expressions/include/GradExpression.h index 58a6c48..cadc85c 100644 --- a/libs/expressions/include/GradExpression.h +++ b/libs/expressions/include/GradExpression.h @@ -39,7 +39,7 @@ namespace cfem::sym std::shared_ptr m_expr; std::shared_ptr createGradient() const override{ - CFEM_ERROR << "GradExpression does not yet support grad(). Cannot compute spatial derivatives for GradExpression."; + CFEM_THROW("GradExpression does not yet support grad(). Cannot compute spatial derivatives for GradExpression."); return nullptr; } std::shared_ptr createTimeDerivative() const override{ @@ -50,7 +50,7 @@ namespace cfem::sym explicit GradExpression(std::shared_ptr expr) : m_expr(std::move(expr)) { if (auto lambda = std::dynamic_pointer_cast(expr) ){ if (!(lambda->grad())) - CFEM_ERROR << "LambdaExpressions are currentely not supported by the automatic differentation."; + CFEM_THROW("LambdaExpressions are currentely not supported by the automatic differentation."); } } diff --git a/libs/expressions/include/TensorExpression.h b/libs/expressions/include/TensorExpression.h index 9ea7542..5a19095 100644 --- a/libs/expressions/include/TensorExpression.h +++ b/libs/expressions/include/TensorExpression.h @@ -38,7 +38,7 @@ namespace cfem::sym private: std::shared_ptr m_components[9]; std::shared_ptr createGradient() const override { - CFEM_ERROR << "Gradient is not yet supported for tensor expressions. It would required a 3rd order tensor which is not implemented."; + CFEM_THROW("Gradient is not yet supported for tensor expressions. It would required a 3rd order tensor which is not implemented."); return nullptr; } std::shared_ptr createTimeDerivative() const override; @@ -60,7 +60,7 @@ namespace cfem::sym for (int i = 0; i < 9; ++i) { if (!m_components[i] || m_components[i]->getDimension() != 1) { - CFEM_ERROR << "TensorExpression component " << i << " must be a valid scalar expression."; + CFEM_THROW("TensorExpression component " + std::to_string(i) + " must be a valid scalar expression."); } } diff --git a/libs/expressions/include/UnaryExpression.tpp b/libs/expressions/include/UnaryExpression.tpp index e26bf5b..9ca88f5 100644 --- a/libs/expressions/include/UnaryExpression.tpp +++ b/libs/expressions/include/UnaryExpression.tpp @@ -31,7 +31,7 @@ namespace cfem::sym // --- Dimension Logic & Validation --- if (m_op == MathOp::Transpose || m_op == MathOp::Trace ) { if (m_dimIn != 9) { - CFEM_ERROR << "Transpose/Trace are only defined for dim 9 tensors."; + CFEM_THROW("Transpose/Trace are only defined for dim 9 tensors."); } } @@ -143,7 +143,7 @@ namespace cfem::sym for(int i=0; igetDimension() != 1) { - CFEM_ERROR << "VectorExpression component " << i << " must be a valid scalar expression. Got dim = " << m_components[i]->getDimension(); + CFEM_THROW("VectorExpression component " << i << " must be a valid scalar expression. Got dim = " << m_components[i]->getDimension() ); } } diff --git a/libs/expressions/src/LambdaExpressionBase.cpp b/libs/expressions/src/LambdaExpressionBase.cpp index 13b86c7..9317ec1 100644 --- a/libs/expressions/src/LambdaExpressionBase.cpp +++ b/libs/expressions/src/LambdaExpressionBase.cpp @@ -30,17 +30,17 @@ namespace cfem::sym m_timeDeriv = std::move(timeDeriv); if (m_gradient){ if (m_dimension == 1 && m_gradient->getDimension() != 3){ - CFEM_ERROR << "The gradient of a scalar lambda expression must be a 3D vector."; + CFEM_THROW("The gradient of a scalar lambda expression must be a 3D vector."); } if (m_dimension == 3 && m_gradient->getDimension() != 9){ - CFEM_ERROR << "The gradient of a vector lambda expression must be a 2nd order tensor (dimension 9)."; + CFEM_THROW("The gradient of a vector lambda expression must be a 2nd order tensor (dimension 9)."); } } if (m_timeDeriv){ if (m_timeDeriv->getDimension() != m_dimension){ - CFEM_ERROR << "The time derivative must have the same dimension (" - << m_dimension << ") as the lambda expression."; + CFEM_THROW( "The time derivative must have the same dimension (" + << m_dimension << ") as the lambda expression."); } } } diff --git a/libs/fem/assembly/CMakeLists.txt b/libs/fem/assembly/CMakeLists.txt index 8f0b450..70a2563 100644 --- a/libs/fem/assembly/CMakeLists.txt +++ b/libs/fem/assembly/CMakeLists.txt @@ -15,4 +15,5 @@ target_include_directories(cfem_assembly target_link_libraries(cfem_assembly PUBLIC cfem_petsc_wrappers cfem_utils -) \ No newline at end of file + cfem_openmp_interface +) diff --git a/libs/fem/fefield/CMakeLists.txt b/libs/fem/fefield/CMakeLists.txt index cb4fb62..4af2806 100644 --- a/libs/fem/fefield/CMakeLists.txt +++ b/libs/fem/fefield/CMakeLists.txt @@ -10,11 +10,12 @@ target_include_directories(cfem_fefield ) target_link_libraries(cfem_fefield - PUBLIC + PUBLIC + cfem_openmp_interface cfem_utils cfem_petsc_wrappers cfem_reference_element cfem_expressions cfem_mesh cfem_fespace -) \ No newline at end of file +) diff --git a/libs/fem/fespace/CMakeLists.txt b/libs/fem/fespace/CMakeLists.txt index 1bd4212..0c6728a 100644 --- a/libs/fem/fespace/CMakeLists.txt +++ b/libs/fem/fespace/CMakeLists.txt @@ -12,9 +12,10 @@ target_include_directories(cfem_fespace target_link_libraries(cfem_fespace PUBLIC + cfem_openmp_interface cfem_utils cfem_reference_element cfem_mesh cfem_expressions cfem_fefield -) \ No newline at end of file +) diff --git a/libs/fem/mesh/CMakeLists.txt b/libs/fem/mesh/CMakeLists.txt index 7c1c6be..2368c51 100644 --- a/libs/fem/mesh/CMakeLists.txt +++ b/libs/fem/mesh/CMakeLists.txt @@ -23,4 +23,5 @@ target_include_directories(cfem_mesh target_link_libraries(cfem_mesh PUBLIC cfem_utils -) \ No newline at end of file + cfem_openmp_interface +) diff --git a/libs/prepostprocessing/numerical_integration/CMakeLists.txt b/libs/prepostprocessing/numerical_integration/CMakeLists.txt index 9827f68..52657b7 100644 --- a/libs/prepostprocessing/numerical_integration/CMakeLists.txt +++ b/libs/prepostprocessing/numerical_integration/CMakeLists.txt @@ -8,9 +8,10 @@ target_include_directories(cfem_numerical_integration target_link_libraries(cfem_numerical_integration PUBLIC cfem_utils + cfem_openmp_interface cfem_reference_element cfem_quadrature cfem_mesh cfem_fespace cfem_expressions -) \ No newline at end of file +) diff --git a/libs/prepostprocessing/numerical_integration/src/FunctionIntegrator.cpp b/libs/prepostprocessing/numerical_integration/src/FunctionIntegrator.cpp index d14a165..78d9311 100644 --- a/libs/prepostprocessing/numerical_integration/src/FunctionIntegrator.cpp +++ b/libs/prepostprocessing/numerical_integration/src/FunctionIntegrator.cpp @@ -35,7 +35,7 @@ namespace cfem m_schemes_params(schemes_params) { if (!m_mesh) - CFEM_ERROR << "FunctionIntegrator requires a valid mesh."; + CFEM_THROW("FunctionIntegrator requires a valid mesh."); } void FunctionIntegrator::L2Norm(const std::shared_ptr &expr, @@ -133,7 +133,7 @@ namespace cfem const MeshEntity& volume_entity) { if( ! m_mesh->isPartOfTheMeshEntities(boundary_entity.getId()) ){ - CFEM_ERROR << boundary_entity.getName() << " is not part of the mesh."; + CFEM_THROW(boundary_entity.getName() << " is not part of the mesh."); } CfemInt entityDim = boundary_entity.getDimension(); @@ -143,23 +143,23 @@ namespace cfem bool isManifoldSurface = (entityDim == 2 && meshDim == 2); if (!isBoundary && !isManifoldSurface) { - CFEM_ERROR << "computeFlux: Flux can only be computed on a boundary (dim " + CFEM_THROW("computeFlux: Flux can only be computed on a boundary (dim " << meshDim - 1 << ") or a surface manifold (dim 2). Target entity '" - << boundary_entity.getName() << "' has dimension " << entityDim; + << boundary_entity.getName() << "' has dimension " << entityDim); } if (isBoundary) { if( ! m_mesh->isPartOfTheMeshEntities(volume_entity.getId()) ){ - CFEM_ERROR << volume_entity.getName() << " is not part of the mesh."; + CFEM_THROW(volume_entity.getName() << " is not part of the mesh."); } if( volume_entity.getDimension() != meshDim) { - CFEM_ERROR << "The body must have the same dimension as the mesh (" << meshDim << ")."; + CFEM_THROW("The body must have the same dimension as the mesh (" << meshDim << ")."); } } if (expr->getDimension() != 3) { - CFEM_ERROR << "computeFlux: The expression must be a 3D vector. " - << "Got dimension " << expr->getDimension(); + CFEM_THROW("computeFlux: The expression must be a 3D vector. " + << "Got dimension " << expr->getDimension()); } auto fluxExpr = sym::dot(expr, sym::geometric_normal()); @@ -211,8 +211,8 @@ namespace cfem if (!targets[i].first) continue; if (targets[i].first->getDimension() != 1) { - CFEM_ERROR << "FunctionIntegrator only accepts scalar expressions. " - << "Please pass each component separately."; + CFEM_THROW( "FunctionIntegrator only accepts scalar expressions. " + << "Please pass each component separately."); } targets[i].first->prepare(); @@ -408,8 +408,8 @@ namespace cfem if (!targets[i].first) continue; if (targets[i].first->getDimension() != 1) { - CFEM_ERROR << "FunctionIntegrator only accepts scalar expressions. " - << "Please pass each component separately."; + CFEM_THROW( "FunctionIntegrator only accepts scalar expressions. " + << "Please pass each component separately."); } targets[i].first->prepare(); diff --git a/libs/utils/include/CfemLogger.h b/libs/utils/include/CfemLogger.h index 40dfecd..abedfd4 100644 --- a/libs/utils/include/CfemLogger.h +++ b/libs/utils/include/CfemLogger.h @@ -20,27 +20,43 @@ #include #include #include -#include +#include +#include +#include "ContractException.h" // Intégration de tes contrats namespace cfem { /** - * @brief Custom exception for CFEM errors + * @brief État global des erreurs (Lock-Free / HPC-Ready) */ -class CfemException : public std::runtime_error { +class ErrorManager { public: - explicit CfemException(const std::string& msg) : std::runtime_error(msg) {} + static ErrorManager& instance() { + static ErrorManager inst; + return inst; + } + + void addWarning() noexcept { warnings_.fetch_add(1, std::memory_order_relaxed); } + void addError() noexcept { errors_.fetch_add(1, std::memory_order_relaxed); } + + bool hasErrors() const noexcept { return errors_.load(std::memory_order_acquire) > 0; } + uint32_t errorCount() const noexcept { return errors_.load(std::memory_order_relaxed); } + + void clear() noexcept { + errors_.store(0); + warnings_.store(0); + } + +private: + ErrorManager() = default; + std::atomic errors_{0}; + std::atomic warnings_{0}; }; -enum class LogLevel { Info, Warning, Error }; +enum class LogLevel { Info, Warning, Error, Fatal }; class Logger { public: - static Logger& getInstance() { - static Logger instance; - return instance; - } - struct LogStream { std::stringstream ss; LogLevel level; @@ -49,64 +65,59 @@ class Logger { LogStream(LogLevel l, const char* f, int ln) : level(l), file(f), line(ln) {} - /** - * @brief Destructor handles printing and throwing exceptions for Error level - */ - ~LogStream() noexcept(false) { // Note: noexcept(false) allows throwing from destructor - std::string prefix; - std::string color; - bool showMetadata = false; -#ifndef NDEBUG - showMetadata = true; -#endif + ~LogStream() { + std::string msg = ss.str(); + const char* color = ""; + const char* prefix = ""; switch (level) { - case LogLevel::Info: - prefix = "[INFO] "; - color = "\033[1;32m"; // Green - showMetadata = false; - break; - case LogLevel::Warning: - prefix = "[WARNING] "; - color = "\033[1;33m"; // Yellow + case LogLevel::Info: + prefix = "[INFO] "; + color = "\033[32m"; break; - case LogLevel::Error: + case LogLevel::Warning: + prefix = "[WARN] "; + color = "\033[33m"; + ErrorManager::instance().addWarning(); break; + case LogLevel::Error: prefix = "[ERROR] "; - color = "\033[1;31m"; // Red - break; + color = "\033[31m"; + ErrorManager::instance().addError(); break; + case LogLevel::Fatal: + prefix = "[FATAL] "; + color = "\033[1;31m"; break; } - - // Output the message - std::cerr << color << prefix << "\033[0m" << ss.str(); - - if (showMetadata) { + + std::cerr << color << prefix << "\033[0m" << msg; +#ifndef NDEBUG + if (level != LogLevel::Info) std::cerr << " (" << file << ":" << line << ")"; - } - +#endif std::cerr << std::endl; - if (level == LogLevel::Error) { - throw CfemException(ss.str()); - } + if (level == LogLevel::Fatal) std::abort(); } - template - LogStream& operator<<(const T& msg) { - ss << msg; - return *this; - } + template + LogStream& operator<<(const T& v) { ss << v; return *this; } }; - LogStream build(LogLevel level, const char* file, int line) { - return LogStream(level, file, line); - } - -private: - Logger() = default; + static LogStream build(LogLevel l, const char* f, int ln) { return LogStream(l, f, ln); } }; } // namespace cfem -#define CFEM_INFO cfem::Logger::getInstance().build(cfem::LogLevel::Info, __FILE__, __LINE__) -#define CFEM_WARN cfem::Logger::getInstance().build(cfem::LogLevel::Warning, __FILE__, __LINE__) -#define CFEM_ERROR cfem::Logger::getInstance().build(cfem::LogLevel::Error, __FILE__, __LINE__) \ No newline at end of file +// --- Macros de Logging (Non-bloquant pour Error/Warning) +#define CFEM_INFO cfem::Logger::build(cfem::LogLevel::Info, __FILE__, __LINE__) +#define CFEM_WARN cfem::Logger::build(cfem::LogLevel::Warning, __FILE__, __LINE__) +#define CFEM_ERROR cfem::Logger::build(cfem::LogLevel::Error, __FILE__, __LINE__) +#define CFEM_FATAL cfem::Logger::build(cfem::LogLevel::Fatal, __FILE__, __LINE__) + +#define CFEM_THROW(msg) \ + do { \ + std::ostringstream ss_; \ + ss_ << msg; \ + cfem::ErrorManager::instance().addError(); \ + std::string tmp_msg = ss_.str(); \ + throw cfem::ContractException(__FILE__, __LINE__, "CFEM_THROW", "Runtime Error", tmp_msg.c_str()); \ + } while(0) diff --git a/libs/utils/include/FixedMatrix.h b/libs/utils/include/FixedMatrix.h index bbf7c16..8b8d15b 100644 --- a/libs/utils/include/FixedMatrix.h +++ b/libs/utils/include/FixedMatrix.h @@ -44,27 +44,27 @@ namespace cfem { * @tparam R Number of rows. * @tparam C Number of columns. */ -template +template struct FixedMatrix { /** * @brief Flat array storage to ensure memory contiguity. */ - T data[R * C] = {0}; + T data[R * C] = {}; /** * @brief Default constructor. Initializes all entries to zero. */ - FixedMatrix() = default; + constexpr FixedMatrix() noexcept = default; /** * @brief Initializer list constructor. * @param list Values provided in row-major order. * @throw CFEM_ASSERT if list size does not match R*C. */ - FixedMatrix(std::initializer_list list) { - CFEM_ASSERT(list.size() == R * C); - int i = 0; - for (auto val : list) data[i++] = val; + constexpr FixedMatrix(std::initializer_list list) { + CFEM_ASSERT(list.size() == R * C, "Invalid initializer_list size."); + CfemInt i = 0; + for (const auto& val : list) data[i++] = val; } /** @@ -73,8 +73,8 @@ struct FixedMatrix { * @param c Column index. * @return Reference to the element at (r, c). */ - inline T& operator()(int r, int c) { - CFEM_ASSERT(r < R && c < C, "Index outof bounds."); + constexpr inline T& operator()(CfemInt r, CfemInt c) { + CFEM_ASSERT(r >= 0 && r < R && c >= 0 && c < C, "Matrix Index out of bounds."); return data[r * C + c]; } @@ -84,41 +84,41 @@ struct FixedMatrix { * @param c Column index. * @return Constant reference to the element at (r, c). */ - inline const T& operator()(int r, int c) const { - CFEM_ASSERT(r < R && c < C); + constexpr inline const T& operator()(CfemInt r, CfemInt c) const { + CFEM_ASSERT(r >= 0 && r < R && c >= 0 && c < C, "Matrix Index out of bounds."); return data[r * C + c]; } /** * @brief Returns number of rows. */ - static constexpr int rows() noexcept { return R; } + [[nodiscard]] static constexpr CfemInt rows() noexcept { return R; } /** * @brief Returns number of columns. */ - static constexpr int cols() noexcept { return C; } + [[nodiscard]] static constexpr CfemInt cols() noexcept { return C; } /** * @brief Returns dimensions as a std::pair {Rows, Cols}. */ - static constexpr std::pair dimensions() noexcept { return {R, C}; } + [[nodiscard]] static constexpr std::pair dimensions() noexcept { return {R, C}; } // --- Basic Operations --- /** @brief Resets all matrix entries to zero. */ - void setZero() { - for (int i = 0; i < R * C; ++i) data[i] = static_cast(0); + constexpr void setZero() noexcept{ + std::fill_n(data, R*C, T{}); } /** * @brief Sets the matrix to Identity (\f$ I \f$). * @note Only available for square matrices. */ - void setIdentity() { - static_assert(R == C, "Identity matrix must be square."); + constexpr void setIdentity() { + CFEM_STATIC_ASSERT(R == C, "Identity matrix must be square."); setZero(); - for (int i = 0; i < R; ++i) (*this)(i, i) = static_cast(1); + for (CfemInt i = 0; i < R; ++i) (*this)(i, i) = static_cast(1); } // --- Arithmetic --- @@ -128,7 +128,7 @@ struct FixedMatrix { */ FixedMatrix operator*(T scalar) const { FixedMatrix res; - for (int i = 0; i < R * C; ++i) res.data[i] = data[i] * scalar; + for (CfemInt i = 0; i < R * C; ++i) res.data[i] = data[i] * scalar; return res; } @@ -138,7 +138,7 @@ struct FixedMatrix { FixedMatrix operator/(T scalar) const { CFEM_ASSERT(std::abs(scalar) > cfem::ZERO_TOLERANCE); FixedMatrix res; - for (int i = 0; i < R * C; ++i) res.data[i] = data[i] / scalar; + for (CfemInt i = 0; i < R * C; ++i) res.data[i] = data[i] / scalar; return res; } @@ -147,7 +147,7 @@ struct FixedMatrix { */ FixedMatrix operator+(const FixedMatrix& other) const { FixedMatrix res; - for (int i = 0; i < R * C; ++i) res.data[i] = data[i] + other.data[i]; + for (CfemInt i = 0; i < R * C; ++i) res.data[i] = data[i] + other.data[i]; return res; } @@ -156,7 +156,7 @@ struct FixedMatrix { */ FixedMatrix operator-(const FixedMatrix& other) const { FixedMatrix res; - for (int i = 0; i < R * C; ++i) res.data[i] = data[i] - other.data[i]; + for (CfemInt i = 0; i < R * C; ++i) res.data[i] = data[i] - other.data[i]; return res; } @@ -164,7 +164,7 @@ struct FixedMatrix { * @brief In-place scalar multiplication. */ FixedMatrix& operator*=(T scalar) { - for (int i = 0; i < R * C; ++i) data[i] *= scalar; + for (CfemInt i = 0; i < R * C; ++i) data[i] *= scalar; return *this; } @@ -173,7 +173,7 @@ struct FixedMatrix { */ FixedMatrix& operator/=(T scalar) { CFEM_ASSERT(std::abs(scalar) > cfem::ZERO_TOLERANCE); - for (int i = 0; i < R * C; ++i) data[i] /= scalar; + for (CfemInt i = 0; i < R * C; ++i) data[i] /= scalar; return *this; } @@ -181,7 +181,7 @@ struct FixedMatrix { * @brief In-place matrix addition. */ FixedMatrix& operator+=(const FixedMatrix& other) { - for (int i = 0; i < R * C; ++i) data[i] += other.data[i]; + for (CfemInt i = 0; i < R * C; ++i) data[i] += other.data[i]; return *this; } @@ -189,7 +189,7 @@ struct FixedMatrix { * @brief In-place matrix subtraction. */ FixedMatrix& operator-=(const FixedMatrix& other) { - for (int i = 0; i < R * C; ++i) data[i] -= other.data[i]; + for (CfemInt i = 0; i < R * C; ++i) data[i] -= other.data[i]; return *this; } @@ -199,8 +199,8 @@ struct FixedMatrix { */ FixedMatrix transpose() const { FixedMatrix res; - for (int i = 0; i < R; ++i) - for (int j = 0; j < C; ++j) + for (CfemInt i = 0; i < R; ++i) + for (CfemInt j = 0; j < C; ++j) res(j, i) = (*this)(i, j); return res; } @@ -210,12 +210,12 @@ struct FixedMatrix { * @brief Matrix-Vector Product: \f$ y = A \cdot x \f$ * @relates FixedMatrix */ -template +template inline FixedVector operator*(const FixedMatrix& m, const FixedVector& v) { FixedVector res; - for (int i = 0; i < R; ++i) { + for (CfemInt i = 0; i < R; ++i) { T sum = 0; - for (int j = 0; j < C; ++j) { + for (CfemInt j = 0; j < C; ++j) { sum += m(i, j) * v[j]; } res[i] = sum; @@ -228,14 +228,14 @@ inline FixedVector operator*(const FixedMatrix& m, const FixedVec * Optimized with ikj loop order for better cache performance. * @relates FixedMatrix */ -template +template inline FixedMatrix operator*(const FixedMatrix& A, const FixedMatrix& B) { FixedMatrix res; res.setZero(); - for (int i = 0; i < R; ++i) { - for (int k = 0; k < Inner; ++k) { + for (CfemInt i = 0; i < R; ++i) { + for (CfemInt k = 0; k < Inner; ++k) { T val = A(i, k); - for (int j = 0; j < C2; ++j) { + for (CfemInt j = 0; j < C2; ++j) { res(i, j) += val * B(k, j); } } @@ -257,7 +257,7 @@ inline Vector3D operator*(const FixedMatrix& A, const Vector3D& } /** @brief Left-side scalar multiplication. @relates FixedMatrix */ -template +template inline FixedMatrix operator*(T scalar, const FixedMatrix& m) { return m * scalar; } @@ -315,11 +315,11 @@ inline FixedMatrix inverse(const FixedMatrix& m, T determinant /** * @brief Prints the matrix in a readable grid format. */ -template +template inline std::ostream& operator<<(std::ostream& os, const FixedMatrix& m) { - for (int i = 0; i < R; ++i) { + for (CfemInt i = 0; i < R; ++i) { os << "| "; - for (int j = 0; j < C; ++j) { + for (CfemInt j = 0; j < C; ++j) { os << m(i, j) << (j == C - 1 ? "" : ", "); } os << " |\n"; diff --git a/libs/utils/include/Vector3D.h b/libs/utils/include/Vector3D.h index 43468e8..d78b7ce 100644 --- a/libs/utils/include/Vector3D.h +++ b/libs/utils/include/Vector3D.h @@ -6,13 +6,13 @@ // --- // --- This software is protected by international copyright laws. // --- Use, distribution, or modification of this software in any form, -// --- source or binary, for personal, academic, or non-commercial -// --- purposes is permitted free of charge, provided that this -// --- copyright notice and this permission notice appear in all +// --- source or binary, for personal, academic, or non-commercial +// --- purposes is permitted free of charge, provided that this +// --- copyright notice and this permission notice appear in all // --- copies and supporting documentation. // --- -// --- The authors and contributors provide this code "as is" without -// --- any express or implied warranty. In no event shall they be +// --- The authors and contributors provide this code "as is" without +// --- any express or implied warranty. In no event shall they be // --- held liable for any damages arising from the use of this software. //************************************************************************ @@ -22,9 +22,9 @@ * @brief Lightweight 3D vector class for Finite Element geometric operations. * @version 0.1 * @date 2026-01-27 - * + * * @copyright Copyright (c) 2026 - * + * */ #ifndef CFEM_VECTOR3D @@ -33,185 +33,369 @@ #include "cfemutils.h" #include #include +#include + +namespace cfem +{ + + /** + * @struct Vector3D + * @brief Represents a vector or point in 3D Euclidean space. + * + * This struct provides the basic linear algebra and vector calculus tools + * required for mesh transformations, Jacobian calculations, and normal extraction. + * It is highly optimized, utilizing `constexpr` for compile-time evaluation where possible. + * + * @par Example Usage: + * \code{.cpp} + * cfem::Vector3D a(1.0, 0.0, 0.0); + * cfem::Vector3D b(0.0, 1.0, 0.0); + * cfem::Vector3D normal = a.cross(b).normalized(); // Returns (0, 0, 1) + * CfemReal area = a.cross(b).norm(); // Returns 1.0 + * \endcode + */ + struct Vector3D + { + CfemReal x = 0.0; ///< X-coordinate + CfemReal y = 0.0; ///< Y-coordinate + CfemReal z = 0.0; ///< Z-coordinate -namespace cfem{ + static constexpr int size = 3; -/** - * @struct Vector3D - * @brief Represents a vector or point in 3D Euclidean space. - * * This struct provides the basic linear algebra and vector calculus tools - * required for mesh transformations, Jacobian calculations, and normal extraction. - */ -struct Vector3D { - CfemReal x = 0.0; ///< X-coordinate - CfemReal y = 0.0; ///< Y-coordinate - CfemReal z = 0.0; ///< Z-coordinate + // --- Basic Arithmetic Operators --- - static constexpr int size = 3; + /** @brief Vector addition. */ + [[nodiscard]] constexpr Vector3D operator+(const Vector3D &other) const noexcept + { + return {x + other.x, y + other.y, z + other.z}; + } - // --- Basic Arithmetic Operators --- + /** + * @brief Vector subtraction. + */ + [[nodiscard]] constexpr Vector3D operator-(const Vector3D &other) const noexcept + { + return {x - other.x, y - other.y, z - other.z}; + } - /** @brief Vector addition. */ - Vector3D operator+(const Vector3D& other) const { - return {x + other.x, y + other.y, z + other.z}; - } + /** @brief Right-side scalar multiplication. */ + [[nodiscard]] constexpr Vector3D operator*(CfemReal scalar) const noexcept + { + return {x * scalar, y * scalar, z * scalar}; + } - /** @brief Vector subtraction. */ - Vector3D operator-(const Vector3D& other) const { - return {x - other.x, y - other.y, z - other.z}; - } + /** + * @brief Right-side scalar division. + * @param scalar The denominator. + * @return A new vector divided by the scalar. + * @note If the scalar is near zero (below cfem::ZERO_TOLERANCE), a CFEM_ERROR is logged, + * but the division proceeds (which may result in Inf or NaN values). + */ + [[nodiscard]] constexpr Vector3D operator/(CfemReal scalar) const + { + if (std::abs(scalar) < cfem::ZERO_TOLERANCE) + CFEM_ERROR << "Division by zero in Vector3D operator/"; + return {x / scalar, y / scalar, z / scalar}; + } - /** @brief Right-side scalar multiplication. */ - Vector3D operator*(CfemReal scalar) const { - return {x * scalar, y * scalar, z * scalar}; - } + // --- Compound Assignment Operators --- - /** @brief Scalar division with zero-tolerance check. */ - Vector3D operator/(CfemReal scalar) const { - if (std::abs(scalar) < cfem::ZERO_TOLERANCE) - CFEM_ERROR << "Division by zero in Vector3D operator/"; - return {x / scalar, y / scalar, z / scalar}; - } + constexpr Vector3D &operator+=(const Vector3D &other) noexcept + { + x += other.x; + y += other.y; + z += other.z; + return *this; + } - // --- Compound Assignment Operators --- + constexpr Vector3D &operator-=(const Vector3D &other) noexcept + { + x -= other.x; + y -= other.y; + z -= other.z; + return *this; + } - Vector3D& operator+=(const Vector3D& other) { - x += other.x; y += other.y; z += other.z; - return *this; - } + /** + * @brief In place multiplication by a scalar + * + * @param scalar + * @return Vector3D& + */ + constexpr Vector3D &operator*=(CfemReal scalar) noexcept + { + x *= scalar; + y *= scalar; + z *= scalar; + return *this; + } - Vector3D& operator-=(const Vector3D& other) { - x -= other.x; y -= other.y; z -= other.z; - return *this; - } + /** + * @brief Evaluates strict bitwise equality. + * @param other The vector to compare with. + * @return True if all components match exactly, false otherwise. + * @warning Do NOT use this for geometric or physical checks due to floating-point + * inaccuracies. Use isApprox() instead for tolerance-based equality. + */ + [[nodiscard]] constexpr bool operator==(const Vector3D& other) const noexcept { + return x == other.x && y == other.y && z == other.z; + } - Vector3D& operator*=(CfemReal scalar) { - x *= scalar; y *= scalar; z *= scalar; - return *this; - } + [[nodiscard]] bool operator!=(const Vector3D &other) const noexcept + { + return !(*this == other); + } - bool operator==(const Vector3D& other) const { - return cfem::NumericPolicy::areEqual(x, other.x) && - cfem::NumericPolicy::areEqual(y, other.y) && - cfem::NumericPolicy::areEqual(z, other.z); - } + /** + * @brief Checks if two vectors are geometrically identical up to a numerical tolerance. + * + * This method relies on the global NumericPolicy to safely compare floating-point + * coordinates. It dynamically computes a local scale for each axis to prevent + * floating-point traps when dealing with very large or very small coordinates. + * + * @param other The vector to compare with. + * @return True if the vectors are identical within the defined tolerance, false otherwise. + */ + [[nodiscard]] bool isApprox(const Vector3D& other) const noexcept { + CfemReal scale_x = std::max(std::abs(x), std::abs(other.x)); + CfemReal scale_y = std::max(std::abs(y), std::abs(other.y)); + CfemReal scale_z = std::max(std::abs(z), std::abs(other.z)); + + return NumericPolicy::areEqual(x, other.x, scale_x) && + NumericPolicy::areEqual(y, other.y, scale_y) && + NumericPolicy::areEqual(z, other.z, scale_z); + } - bool operator!=(const Vector3D& other) const { - return !(*this == other); - } + /** + * @brief In-place Scalar division. + * @param scalar The denominator. + * @note If the scalar is near zero (below cfem::ZERO_TOLERANCE), a CFEM_ERROR is logged, + * but the division proceeds (which may result in Inf or NaN values). + */ + Vector3D &operator/=(CfemReal scalar) + { + if (std::abs(scalar) < cfem::ZERO_TOLERANCE) + CFEM_ERROR << "Division by zero in Vector3D operator/="; + x /= scalar; + y /= scalar; + z /= scalar; + return *this; + } - Vector3D& operator/=(CfemReal scalar) { - if (std::abs(scalar) < cfem::ZERO_TOLERANCE) - CFEM_ERROR << "Division by zero in Vector3D operator/="; - x /= scalar; y /= scalar; z /= scalar; - return *this; - } - - // --- Accessors --- - - /** @brief Read-only access by index (0=x, 1=y, 2=z). */ - CfemReal operator[](int i) const { - if (i == 0) return x; - if (i == 1) return y; - if (i == 2) return z; - CFEM_ERROR << "Index out of bound in Vector3D::operator[] const."; - return 0.0; - } + // --- Accessors --- + + /** @brief Read-only access by index (0=x, 1=y, 2=z). */ + [[nodiscard]] CfemReal operator[](int i) const + { + if (i == 0) + return x; + if (i == 1) + return y; + if (i == 2) + return z; + CFEM_THROW("Index out of bound in Vector3D::operator[] const."); + return 0.0; + } - /** @brief Read-write access by index (0=x, 1=y, 2=z). */ - CfemReal& operator[](int i) { - if (i == 0) return x; - if (i == 1) return y; - if (i == 2) return z; - CFEM_ERROR << "Index out of bound in Vector3D::operator[]."; - return x; - } - - // --- Vector Calculus --- + /** + * @brief Read-write access by index (0=x, 1=y, 2=z). + */ + CfemReal &operator[](int i) + { + if (i == 0) + return x; + if (i == 1) + return y; + if (i == 2) + return z; + CFEM_THROW("Index out of bound in Vector3D::operator[]."); + return x; + } - /** * @brief Computes the dot (inner) product. - * @return a.x*b.x + a.y*b.y + a.z*b.z - */ - CfemReal dot(const Vector3D& other) const { - return x * other.x + y * other.y + z * other.z; - } + // --- Vector Calculus --- - /** * @brief Computes the cross product. - * @return A vector perpendicular to both this and other. - */ - Vector3D cross(const Vector3D& other) const { - return { - y * other.z - z * other.y, - z * other.x - x * other.z, - x * other.y - y * other.x - }; - } - - /** @brief Returns the squared Euclidean norm (faster than norm()). */ - CfemReal normSquare() const { - return x * x + y * y + z * z; + /** * @brief Computes the dot (inner) product. + * @return a.x*b.x + a.y*b.y + a.z*b.z + */ + [[nodiscard]] constexpr CfemReal dot(const Vector3D &other) const noexcept + { + return x * other.x + y * other.y + z * other.z; + } + + /** * @brief Computes the cross product. + * @return A vector perpendicular to both this and other. + */ + [[nodiscard]] constexpr Vector3D cross(const Vector3D &other) const noexcept + { + return { + y * other.z - z * other.y, + z * other.x - x * other.z, + x * other.y - y * other.x}; + } + + /** + * @brief Computes the scalar triple product: this . (b x c). + * Extremely useful for computing 3D element volumes and Jacobian determinants. + */ + [[nodiscard]] constexpr CfemReal scalarTripleProduct(const Vector3D& b, const Vector3D& c) const { + return this->dot(b.cross(c)); + } + + /** @brief Returns the squared Euclidean norm (faster than norm()). */ + [[nodiscard]] constexpr inline CfemReal normSquare() const noexcept + { + return x * x + y * y + z * z; + } + + /** @brief Returns the Euclidean norm (length) of the vector. */ + [[nodiscard]] inline CfemReal norm() const noexcept + { + return std::sqrt(normSquare()); + } + + /** + * @brief Returns a unit vector in the same direction. + * Returns {0,0,0} if the vector is near-null. + */ + [[nodiscard]] inline Vector3D normalized() const + { + CfemReal n = norm(); + if (n < cfem::ZERO_TOLERANCE) + return {0, 0, 0}; + return (*this) / n; + } + + inline Vector3D &normalize() noexcept + { + CfemReal n = norm(); + if (n > cfem::ZERO_TOLERANCE) + { + x /= n; + y /= n; + z /= n; + } + return *this; + } + + /** + * @brief Returns the squared distance to another point. (Faster, no sqrt) + */ + [[nodiscard]] constexpr CfemReal distanceSquaredTo(const Vector3D& other) const { + return (*this - other).normSquare(); + } + + /** + * @brief Returns the Euclidean distance to another point. + */ + [[nodiscard]] CfemReal distanceTo(const Vector3D& other) const { + return std::sqrt(distanceSquaredTo(other)); + } + + [[nodiscard]] inline bool isClose(const Vector3D &other, CfemReal tol = cfem::ZERO_TOLERANCE) const noexcept + { + return (std::abs(x - other.x) < tol && + std::abs(y - other.y) < tol && + std::abs(z - other.z) < tol); + } + + /** + * @brief Computes the angle (in radians) between this vector and another. + * Safe against floating point inaccuracies. + */ + [[nodiscard]] inline CfemReal angleWith(const Vector3D& other) const { + return std::acos(cosAngleWith(other)); + } + + /** + * @brief Computes the cosine of the angle between this vector and another. + * @param other The second vector. + * @return The cosine of the angle in the range [-1.0, 1.0]. + * @note Faster than angleWith() and numerically stable. + */ + [[nodiscard]] CfemReal cosAngleWith(const Vector3D& other) const { + CfemReal denom = this->norm() * other.norm(); + + if (denom < cfem::ZERO_TOLERANCE) return 0.0; + + return std::clamp(this->dot(other) / denom, -1.0, 1.0); // std::clamp to avoid floating point errors (the result should stay between -1 and 1) + } + + /** + * @brief Computes the orthogonal projection of this vector onto a target direction. + * @param target The vector defining the projection direction. + * @return The projected Vector3D. Returns {0,0,0} if the target is a null vector. + */ + [[nodiscard]] Vector3D projectOnto(const Vector3D& target) const { + CfemReal targetNormSq = target.normSquare(); + if (targetNormSq < cfem::ZERO_TOLERANCE) return {0, 0, 0}; + return target * (this->dot(target) / targetNormSq); + } + + /** + * @brief Returns the component of this vector perpendicular to 'target'. + */ + [[nodiscard]] Vector3D perpendicularTo(const Vector3D& target) const { + return *this - projectOnto(target); + } + + /** + * @brief Returns a vector containing the minimum components of both vectors. + */ + [[nodiscard]] constexpr Vector3D cwiseMin(const Vector3D& other) const { + return {std::min(x, other.x), std::min(y, other.y), std::min(z, other.z)}; + } + + /** + * @brief Returns a vector containing the maximum components of both vectors. + */ + [[nodiscard]] constexpr Vector3D cwiseMax(const Vector3D& other) const { + return {std::max(x, other.x), std::max(y, other.y), std::max(z, other.z)}; + } + + /** + * @brief Returns the vector resulting from component-wise multiplication. + */ + [[nodiscard]] constexpr Vector3D cwiseProduct(const Vector3D& other) const { + return {x * other.x, y * other.y, z * other.z}; + } + }; + + // --- Non-member Helper Functions --- + + /** @brief Left-side scalar multiplication (e.g., 2.0 * vector). */ + [[nodiscard]] inline constexpr Vector3D operator*(CfemReal scalar, const Vector3D &p) noexcept + { + return {p.x * scalar, p.y * scalar, p.z * scalar}; } /** @brief Returns the Euclidean norm (length) of the vector. */ - CfemReal norm() const { - return std::sqrt(normSquare()); + [[nodiscard]] inline CfemReal norm(const Vector3D &v) noexcept + { + return v.norm(); } - /** * @brief Returns a unit vector in the same direction. - * Returns {0,0,0} if the vector is near-null. - */ - Vector3D normalized() const { - CfemReal n = norm(); - if (n < cfem::ZERO_TOLERANCE) return {0, 0, 0}; - return (*this) / n; + /** @brief Functional style dot product. */ + [[nodiscard]] inline constexpr CfemReal dot(const Vector3D &a, const Vector3D &b) noexcept + { + return a.dot(b); } - const Vector3D& normalize() { - CfemReal n = norm(); - if (n > cfem::ZERO_TOLERANCE) { - x /= n; y /= n; z /= n; - } - return *this; + /** @brief Functional style cross product. */ + [[nodiscard]] inline constexpr Vector3D cross(const Vector3D &a, const Vector3D &b) noexcept + { + return a.cross(b); } - bool isClose(const Vector3D& other, CfemReal tol = cfem::ZERO_TOLERANCE) const { - return (std::abs(x - other.x) < tol && - std::abs(y - other.y) < tol && - std::abs(z - other.z) < tol); + /** @brief Overloads the output stream operator for easy printing/logging. + * @param os Output stream (e.g., std::cout or a file stream). + * @param p The Vector3D to print. + * @return Reference to the stream for chaining. + */ + inline std::ostream &operator<<(std::ostream &os, const Vector3D &p) + { + return os << "(" << p.x << ", " << p.y << ", " << p.z << ")"; } -}; - -// --- Non-member Helper Functions --- - -/** @brief Left-side scalar multiplication (e.g., 2.0 * vector). */ -inline Vector3D operator*(CfemReal scalar, const Vector3D& p) { - return {p.x * scalar, p.y * scalar, p.z * scalar}; -} - -/** @brief Returns the Euclidean norm (length) of the vector. */ -inline CfemReal norm(const Vector3D& v) { - return v.norm(); -} - -/** @brief Functional style dot product. */ -inline CfemReal dot(const Vector3D& a, const Vector3D& b) { - return a.dot(b); -} - -/** @brief Functional style cross product. */ -inline Vector3D cross(const Vector3D& a, const Vector3D& b) { - return a.cross(b); -} - -/** @brief Overloads the output stream operator for easy printing/logging. - * @param os Output stream (e.g., std::cout or a file stream). - * @param p The Vector3D to print. - * @return Reference to the stream for chaining. - */ -inline std::ostream& operator<<(std::ostream& os, const Vector3D& p) { - return os << "(" << p.x << ", " << p.y << ", " << p.z << ")"; -} -} // namespace +} // namespace #endif // CFEM_VECTOR3D \ No newline at end of file diff --git a/libs/utils/include/cfemutils.h b/libs/utils/include/cfemutils.h index dc2d670..8affdcc 100644 --- a/libs/utils/include/cfemutils.h +++ b/libs/utils/include/cfemutils.h @@ -48,6 +48,7 @@ namespace cfem * is effectively zero. */ constexpr CfemReal ZERO_TOLERANCE = 1e-14; + constexpr CfemReal ZERO_TOLERANCE_SQ = 1e-28; inline constexpr CfemReal pi = static_cast(3.141592653589793238462643383279502884L); /** diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 2652ee3..b6414e7 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -21,6 +21,7 @@ add_executable(unit_tests target_link_libraries(unit_tests PRIVATE GTest::gtest + cfem_openmp_interface cfem_solvers cfem_utils cfem_quadrature @@ -31,4 +32,6 @@ target_link_libraries(unit_tests cfem_numerical_integration ) -gtest_discover_tests(unit_tests) \ No newline at end of file +target_compile_definitions(unit_tests PRIVATE CFEM_TESTING) + +gtest_discover_tests(unit_tests)