diff --git a/kratos/includes/kratos_application.h b/kratos/includes/kratos_application.h index f3ffa3c50243..5bad37c27933 100644 --- a/kratos/includes/kratos_application.h +++ b/kratos/includes/kratos_application.h @@ -81,6 +81,7 @@ #include "modeler/combine_model_part_modeler.h" #include "modeler/connectivity_preserve_modeler.h" #include "modeler/voxel_mesh_generator_modeler.h" +#include "modeler/cartesian_mesh_generator_modeler.h" #include "modeler/clean_up_problematic_triangles_modeler.h" namespace Kratos { @@ -571,6 +572,7 @@ class KRATOS_API(KRATOS_CORE) KratosApplication { const CombineModelPartModeler mCombineModelPartModeler; const ConnectivityPreserveModeler mConnectivityPreserveModeler; const VoxelMeshGeneratorModeler mVoxelMeshGeneratorModeler; + const CartesianMeshGeneratorModeler mCartesianMeshGeneratorModeler; const CleanUpProblematicTrianglesModeler mCleanUpProblematicTrianglesModeler; // Base constitutive law definition diff --git a/kratos/modeler/cartesian_mesh_generator_modeler.cpp b/kratos/modeler/cartesian_mesh_generator_modeler.cpp new file mode 100644 index 000000000000..ba44153808d2 --- /dev/null +++ b/kratos/modeler/cartesian_mesh_generator_modeler.cpp @@ -0,0 +1,541 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Jordi Cotela Dalmau +// + +// System includes + +// External includes + +// Project includes +#include "utilities/timer.h" +#include "modeler/cartesian_mesh_generator_modeler.h" + +namespace Kratos +{ + +CartesianMeshGeneratorModeler::CartesianMeshGeneratorModeler() + : Modeler() + , mpModel(nullptr) + , mpSourceModelPart(nullptr) + , mElementSize(0.0) +{ +} + +/***********************************************************************************/ +/***********************************************************************************/ + +CartesianMeshGeneratorModeler::CartesianMeshGeneratorModeler( + Model& rModel, + Parameters ModelerParameters + ) + : Modeler(rModel, ModelerParameters) + , mpModel(&rModel) + , mpSourceModelPart(nullptr) +{ + mParameters.ValidateAndAssignDefaults(GetDefaultParameters()); + mElementSize = mParameters["element_size"].GetDouble(); +} + +/***********************************************************************************/ +/***********************************************************************************/ + +const Parameters CartesianMeshGeneratorModeler::GetDefaultParameters() const +{ + return Parameters(R"({ + "input_model_part_name" : "", + "output_model_part_name" : "", + "element_name" : "Element3D4N", + "element_size" : 1.0 + })"); +} + +/***********************************************************************************/ +/***********************************************************************************/ + +void CartesianMeshGeneratorModeler::SetupModelPart() +{ + KRATOS_TRY + + KRATOS_ERROR_IF(mpModel == nullptr) + << "CartesianMeshGeneratorModeler::SetupModelPart called without a Model." << std::endl; + + const std::string input_name = mParameters["input_model_part_name"].GetString(); + const std::string output_name = mParameters["output_model_part_name"].GetString(); + const std::string element_name = mParameters["element_name"].GetString(); + + KRATOS_ERROR_IF(input_name.empty()) + << "CartesianMeshGeneratorModeler: \"input_model_part_name\" must not be empty." << std::endl; + KRATOS_ERROR_IF(output_name.empty()) + << "CartesianMeshGeneratorModeler: \"output_model_part_name\" must not be empty." << std::endl; + + mpSourceModelPart = &mpModel->GetModelPart(input_name); + + ModelPart& r_output = mpModel->HasModelPart(output_name) + ? mpModel->GetModelPart(output_name) + : mpModel->CreateModelPart(output_name); + + if (!r_output.RecursivelyHasProperties(0)) + r_output.CreateNewProperties(0); + + const Element& r_ref_element = KratosComponents::Get(element_name); + GenerateCartesianMesh(r_output, r_ref_element); + + KRATOS_CATCH("") +} + +/***********************************************************************************/ +/***********************************************************************************/ + +void CartesianMeshGeneratorModeler::GenerateCartesianMesh( + ModelPart& rModelPart, + Element const& rReferenceElement + ) +{ + KRATOS_TRY + + KRATOS_ERROR_IF(mpSourceModelPart == nullptr) + << "CartesianMeshGeneratorModeler: source model part is not set." << std::endl; + + const unsigned int dimension = rReferenceElement.GetGeometry().WorkingSpaceDimension(); + + KRATOS_INFO("CartesianMeshGeneratorModeler") << "Generating mesh for dimension: " << dimension << std::endl; + + Timer::Start("Generating Mesh"); + + CalculateBoundingBox(*mpSourceModelPart, mMinPoint, mMaxPoint); + CalculateDivisionNumbers(); + + unsigned int start_node_id = mpSourceModelPart->NumberOfNodes() + 1; + unsigned int start_element_id = 1; + + if (dimension == 3) { + // 3D Cartesian tetrahedral mesh using Freudenthal decomposition (6 tets per hex cell) + Timer::Start("Generating Nodes"); + + const double x0 = mMinPoint.X(); + const double y0 = mMinPoint.Y(); + const double z0 = mMinPoint.Z(); + const unsigned int nx = mDivisionsNumber[0]; + const unsigned int ny = mDivisionsNumber[1]; + const unsigned int nz = mDivisionsNumber[2]; + + // 3D index with correct strides (nx, nx*ny) + auto NodeIdx3D = [&](unsigned int i, unsigned int j, unsigned int k) -> unsigned int { + return i + nx * j + nx * ny * k; + }; + + const unsigned int total_nodes = nx * ny * nz; + std::vector node_ptrs(total_nodes); + ModelPart::NodesContainerType::ContainerType& r_nodes_array = rModelPart.NodesArray(); + r_nodes_array.resize(total_nodes); + + auto p_variables_list = rModelPart.pGetNodalSolutionStepVariablesList(); + const SizeType buffer_size = rModelPart.GetBufferSize(); + for (unsigned int k = 0; k < nz; k++) { + for (unsigned int j = 0; j < ny; j++) { + for (unsigned int i = 0; i < nx; i++) { + auto p_node = Node::Pointer(new Node( + start_node_id++, + x0 + i * mElementSize, + y0 + j * mElementSize, + z0 + k * mElementSize)); + p_node->SetSolutionStepVariablesList(p_variables_list); + p_node->SetBufferSize(buffer_size); + node_ptrs[NodeIdx3D(i,j,k)] = p_node; + r_nodes_array[NodeIdx3D(i,j,k)] = p_node; + } + } + } + + Timer::Stop("Generating Nodes"); + Timer::Start("Generating Elements"); + + const unsigned int sx = mSegmentsNumber[0]; + const unsigned int sy = mSegmentsNumber[1]; + const unsigned int sz = mSegmentsNumber[2]; + ModelPart::ElementsContainerType::ContainerType& r_elements_array = rModelPart.ElementsArray(); + r_elements_array.resize(6 * sx * sy * sz); + + Element::NodesArrayType tet_nodes(4); + unsigned int elem_counter = 0; + + for (unsigned int k = 0; k < sz; k++) { + for (unsigned int j = 0; j < sy; j++) { + for (unsigned int i = 0; i < sx; i++) { + auto c000 = node_ptrs[NodeIdx3D(i, j, k )]; + auto c100 = node_ptrs[NodeIdx3D(i+1, j, k )]; + auto c010 = node_ptrs[NodeIdx3D(i, j+1, k )]; + auto c110 = node_ptrs[NodeIdx3D(i+1, j+1, k )]; + auto c001 = node_ptrs[NodeIdx3D(i, j, k+1)]; + auto c101 = node_ptrs[NodeIdx3D(i+1, j, k+1)]; + auto c011 = node_ptrs[NodeIdx3D(i, j+1, k+1)]; + auto c111 = node_ptrs[NodeIdx3D(i+1, j+1, k+1)]; + + // Freudenthal decomposition: 6 tets all sharing the c000-c111 body diagonal + tet_nodes(0)=c000; tet_nodes(1)=c100; tet_nodes(2)=c110; tet_nodes(3)=c111; + r_elements_array[elem_counter++] = rReferenceElement.Create(start_element_id++, tet_nodes, rReferenceElement.pGetProperties()); + + tet_nodes(0)=c000; tet_nodes(1)=c100; tet_nodes(2)=c101; tet_nodes(3)=c111; + r_elements_array[elem_counter++] = rReferenceElement.Create(start_element_id++, tet_nodes, rReferenceElement.pGetProperties()); + + tet_nodes(0)=c000; tet_nodes(1)=c010; tet_nodes(2)=c110; tet_nodes(3)=c111; + r_elements_array[elem_counter++] = rReferenceElement.Create(start_element_id++, tet_nodes, rReferenceElement.pGetProperties()); + + tet_nodes(0)=c000; tet_nodes(1)=c010; tet_nodes(2)=c011; tet_nodes(3)=c111; + r_elements_array[elem_counter++] = rReferenceElement.Create(start_element_id++, tet_nodes, rReferenceElement.pGetProperties()); + + tet_nodes(0)=c000; tet_nodes(1)=c001; tet_nodes(2)=c101; tet_nodes(3)=c111; + r_elements_array[elem_counter++] = rReferenceElement.Create(start_element_id++, tet_nodes, rReferenceElement.pGetProperties()); + + tet_nodes(0)=c000; tet_nodes(1)=c001; tet_nodes(2)=c011; tet_nodes(3)=c111; + r_elements_array[elem_counter++] = rReferenceElement.Create(start_element_id++, tet_nodes, rReferenceElement.pGetProperties()); + } + } + } + + Timer::Stop("Generating Elements"); + Timer::Stop("Generating Mesh"); + return; + } + + KRATOS_ERROR_IF(dimension != 2) << "Only 2D and 3D meshes are supported!" << std::endl; + + // 2D path: quadrilateral elements with inside-test via ray casting + const unsigned int number_of_nodes = mDivisionsNumber[0] * mDivisionsNumber[1] * mDivisionsNumber[2]; + const unsigned int number_of_elements = mSegmentsNumber[0] * mSegmentsNumber[1] * mSegmentsNumber[2]; + + KRATOS_INFO("CartesianMeshGeneratorModeler") << "Number of nodes: " << number_of_nodes << std::endl; + KRATOS_INFO("CartesianMeshGeneratorModeler") << "Number of elements: " << number_of_elements << std::endl; + + CalculateBoundaryIntersections(*mpSourceModelPart); + CalculateIsInside(*mpSourceModelPart); + CalculateNormals(); + + Timer::Start("Generating Nodes"); + + const double x0 = mMinPoint.X(); + const double y0 = mMinPoint.Y(); + const double z0 = mMinPoint.Z(); + + ModelPart::NodesContainerType::ContainerType& r_nodes_array = rModelPart.NodesArray(); + ModelPart::NodesContainerType::ContainerType temp_nodes_array(number_of_nodes); + + ModelPart::ElementsContainerType::ContainerType& r_elements_array = rModelPart.ElementsArray(); + + auto p_variables_list = rModelPart.pGetNodalSolutionStepVariablesList(); + const SizeType buffer_size = rModelPart.GetBufferSize(); + for (unsigned int k = 0; k < mDivisionsNumber[2]; k++) + for (unsigned int j = 0; j < mDivisionsNumber[1]; j++) + for (unsigned int i = 0; i < mDivisionsNumber[0]; i++) { + auto p_node = Node::Pointer(new Node( + start_node_id++, + x0 + i * mElementSize, + y0 + j * mElementSize, + z0 + k * mElementSize)); + p_node->SetSolutionStepVariablesList(p_variables_list); + p_node->SetBufferSize(buffer_size); + temp_nodes_array[NodeIndex(i,j,k)] = p_node; + } + + unsigned int number_of_active_nodes = 0; + for (unsigned int i = 0; i < number_of_nodes; i++) + if (mIsInside[i]) + number_of_active_nodes++; + + r_nodes_array.resize(number_of_active_nodes); + + unsigned int index = 0; + for (unsigned int i = 0; i < number_of_nodes; i++) + if (mIsInside[i]) + r_nodes_array[index++] = temp_nodes_array[i]; + + Timer::Stop("Generating Nodes"); + Timer::Start("Generating Elements"); + + Element::NodesArrayType element_nodes(4); + + unsigned int number_of_active_elements = 0; + for (unsigned int j = 0; j < mSegmentsNumber[1]; j++) + for (unsigned int i = 0; i < mSegmentsNumber[0]; i++) + if (mIsInside[NodeIndex(i,j,0)] & mIsInside[NodeIndex(i+1,j,0)] & + mIsInside[NodeIndex(i+1,j+1,0)] & mIsInside[NodeIndex(i,j+1,0)]) + number_of_active_elements++; + + r_elements_array.resize(number_of_active_elements); + + unsigned int counter = 0; + for (unsigned int j = 0; j < mSegmentsNumber[1]; j++) { + for (unsigned int i = 0; i < mSegmentsNumber[0]; i++) { + if (mIsInside[NodeIndex(i,j,0)] & mIsInside[NodeIndex(i+1,j,0)] & + mIsInside[NodeIndex(i+1,j+1,0)] & mIsInside[NodeIndex(i,j+1,0)]) { + element_nodes(0) = temp_nodes_array[NodeIndex(i, j, 0)]; + element_nodes(1) = temp_nodes_array[NodeIndex(i+1, j, 0)]; + element_nodes(2) = temp_nodes_array[NodeIndex(i+1, j+1, 0)]; + element_nodes(3) = temp_nodes_array[NodeIndex(i, j+1, 0)]; + r_elements_array[counter++] = rReferenceElement.Create(start_element_id++, element_nodes, rReferenceElement.pGetProperties()); + } + } + } + + for (auto& r_elem : mpSourceModelPart->Elements()) { + auto& r_geometry = r_elem.GetGeometry(); + array_1d coords1, coords2; + noalias(coords1) = r_geometry[0].Coordinates() + mNormals[r_geometry[0].Id()-1]; + noalias(coords2) = r_geometry[1].Coordinates() + mNormals[r_geometry[1].Id()-1]; + Point point1(coords1[0], coords1[1], coords1[2]); + Point point2(coords2[0], coords2[1], coords2[2]); + unsigned int index1 = FindNearestNodeIndex(point1, mNormals[r_geometry[0].Id()-1]); + unsigned int index2 = FindNearestNodeIndex(point2, mNormals[r_geometry[1].Id()-1]); + element_nodes(0) = r_geometry(0); + element_nodes(1) = r_geometry(1); + element_nodes(2) = temp_nodes_array[index2]; + element_nodes(3) = temp_nodes_array[index1]; + r_elements_array.push_back(rReferenceElement.Create(start_element_id++, element_nodes, rReferenceElement.pGetProperties())); + } + + Timer::Stop("Generating Elements"); + Timer::Stop("Generating Mesh"); + + KRATOS_CATCH("") +} + +/***********************************************************************************/ +/***********************************************************************************/ + +void CartesianMeshGeneratorModeler::CalculateNormals() +{ + const array_1d zero = ZeroVector(3); + + if (mNormals.size() != mpSourceModelPart->NumberOfNodes()) + mNormals.resize(mpSourceModelPart->NumberOfNodes(), zero); + else + std::fill(mNormals.begin(), mNormals.end(), zero); + + const double coefficient = mElementSize / 2.0; + + for (auto& r_elem : mpSourceModelPart->Elements()) { + auto& r_geometry = r_elem.GetGeometry(); + array_1d normal; + normal[0] = r_geometry[1].Y() - r_geometry[0].Y(); + normal[1] = -(r_geometry[1].X() - r_geometry[0].X()); + normal[2] = 0.0; + normal *= coefficient / r_geometry.Length(); + mNormals[r_geometry[0].Id()-1] += normal; + mNormals[r_geometry[1].Id()-1] += normal; + } + + for (auto& r_node : mpSourceModelPart->Nodes()) + noalias(r_node.FastGetSolutionStepValue(NORMAL)) = mNormals[r_node.Id()-1]; +} + +/***********************************************************************************/ +/***********************************************************************************/ + +unsigned int CartesianMeshGeneratorModeler::FindNearestNodeIndex( + Point& rThisPoint, + array_1d& rNormal + ) +{ + const double x = (rThisPoint.X() - mMinPoint.X()) / mElementSize; + const double y = (rThisPoint.Y() - mMinPoint.Y()) / mElementSize; + + unsigned int i = static_cast(x); + unsigned int j = static_cast(y); + + if (mIsInside[NodeIndex(i,j,0)]) + return NodeIndex(i,j,0); + + if (rNormal[0] >= 0.0) { + if (rNormal[1] >= 0.0) { + if (rNormal[0] > rNormal[1]) i++; + else j++; + } else { + if (rNormal[0] > -rNormal[1]) i++; + else j--; + } + } else { + if (rNormal[1] >= 0.0) { + if (-rNormal[0] > rNormal[1]) i--; + else j++; + } else { + if (-rNormal[0] > -rNormal[1]) i--; + else j--; + } + } + + return NodeIndex(i,j,0); +} + +/***********************************************************************************/ +/***********************************************************************************/ + +void CartesianMeshGeneratorModeler::CalculateIsInside(ModelPart& rModelPart) +{ + const unsigned int number_of_nodes = mDivisionsNumber[0] * mDivisionsNumber[1] * mDivisionsNumber[2]; + + if (mIsInside.size() != number_of_nodes) + mIsInside.resize(number_of_nodes, 0); + else + std::fill(mIsInside.begin(), mIsInside.end(), 0); + + const int size = static_cast(mSegmentsNumber[1]) + 1; + + for (int j = 0; j < size; j++) { + std::vector& r_j_intersections = mIntersections[j]; + + for (auto j_x = r_j_intersections.begin(); j_x != r_j_intersections.end(); j_x++) { + const auto start = j_x++; + const unsigned int i_start = static_cast(*start / mElementSize); + const unsigned int i_end = static_cast(*j_x / mElementSize); + for (unsigned int i = i_start + 1; i < i_end; i++) + mIsInside[NodeIndex(i, j, 0)] = true; + } + } +} + +/***********************************************************************************/ +/***********************************************************************************/ + +void CartesianMeshGeneratorModeler::CalculateBoundaryIntersections(ModelPart& rModelPart) +{ + if (mIntersections.size() != mSegmentsNumber[1] + 1) + mIntersections.resize(mSegmentsNumber[1] + 1); + + for (auto& r_elem : rModelPart.Elements()) { + const auto& r_geometry = r_elem.GetGeometry(); + + double x1, x2, y1, y2; + if (r_geometry[0].Y() < r_geometry[1].Y()) { + x1 = r_geometry[0].X() - mMinPoint.X(); + x2 = r_geometry[1].X() - mMinPoint.X(); + y1 = r_geometry[0].Y() - mMinPoint.Y(); + y2 = r_geometry[1].Y() - mMinPoint.Y(); + } else { + x1 = r_geometry[1].X() - mMinPoint.X(); + x2 = r_geometry[0].X() - mMinPoint.X(); + y1 = r_geometry[1].Y() - mMinPoint.Y(); + y2 = r_geometry[0].Y() - mMinPoint.Y(); + } + + unsigned int i_start = static_cast(y1 / mElementSize); + const unsigned int i_end = static_cast(y2 / mElementSize); + + if (i_start * mElementSize < y1) + i_start++; + + const double m = (y1 != y2) ? (x2 - x1) / (y2 - y1) : 0.0; + const double delta_x = m * mElementSize; + double x = x1 + (i_start * mElementSize - y1) * m; + + for (unsigned int i = i_start; i <= i_end; i++) { + mIntersections[i].push_back(x); + x += delta_x; + } + } + + for (auto& r_row : mIntersections) + std::sort(r_row.begin(), r_row.end()); +} + +/***********************************************************************************/ +/***********************************************************************************/ + +void CartesianMeshGeneratorModeler::CalculateBoundingBox( + ModelPart& rModelPart, + Point& rMinPoint, + Point& rMaxPoint + ) +{ + if (rModelPart.NumberOfElements() == 0 && rModelPart.NumberOfNodes() == 0) { + rMinPoint = Point(); + rMaxPoint = Point(); + return; + } + + // Fall back to node-based BB when no elements are present (e.g. condition-only model parts from STL) + if (rModelPart.NumberOfElements() == 0) { + rMinPoint = *rModelPart.NodesBegin(); + rMaxPoint = rMinPoint; + for (const auto& r_node : rModelPart.Nodes()) { + for (unsigned int i = 0; i < 3; i++) { + if (rMinPoint[i] > r_node[i]) rMinPoint[i] = r_node[i]; + if (rMaxPoint[i] < r_node[i]) rMaxPoint[i] = r_node[i]; + } + } + return; + } + + if (rModelPart.ElementsBegin()->GetGeometry().empty()) { + rMinPoint = Point(); + rMaxPoint = Point(); + return; + } + + rMinPoint = rModelPart.ElementsBegin()->GetGeometry()[0]; + rMaxPoint = rMinPoint; + + for (const auto& r_elem : rModelPart.Elements()) { + for (const auto& r_point : r_elem.GetGeometry()) { + for (unsigned int i = 0; i < Point::Dimension(); i++) { + if (rMinPoint[i] > r_point[i]) rMinPoint[i] = r_point[i]; + if (rMaxPoint[i] < r_point[i]) rMaxPoint[i] = r_point[i]; + } + } + } +} + +/***********************************************************************************/ +/***********************************************************************************/ + +void CartesianMeshGeneratorModeler::CalculateDivisionNumbers() +{ + if (mElementSize == 0.0) + return; + + for (unsigned int i = 0; i < Point::Dimension(); i++) { + const double delta = mMaxPoint[i] - mMinPoint[i]; + int segments_number = static_cast(delta / mElementSize); + + if ((segments_number * mElementSize) < delta) + segments_number++; + + mSegmentsNumber[i] = segments_number; + mDivisionsNumber[i] = segments_number + 1; + + if (mSegmentsNumber[i] == 0) + mSegmentsNumber[i]++; + } +} + +/***********************************************************************************/ +/***********************************************************************************/ + +std::string CartesianMeshGeneratorModeler::Info() const +{ + return "CartesianMeshGeneratorModeler"; +} + +/***********************************************************************************/ +/***********************************************************************************/ + +void CartesianMeshGeneratorModeler::PrintInfo(std::ostream& rOStream) const +{ + rOStream << Info(); +} + +/***********************************************************************************/ +/***********************************************************************************/ + +void CartesianMeshGeneratorModeler::PrintData(std::ostream& rOStream) const +{ +} + +} // namespace Kratos \ No newline at end of file diff --git a/kratos/modeler/cartesian_mesh_generator_modeler.h b/kratos/modeler/cartesian_mesh_generator_modeler.h index 55a9b48cbb71..3ce038d46e78 100644 --- a/kratos/modeler/cartesian_mesh_generator_modeler.h +++ b/kratos/modeler/cartesian_mesh_generator_modeler.h @@ -2,54 +2,54 @@ // ' / __| _` | __| _ \ __| // . \ | ( | | ( |\__ ` // _|\_\_| \__,_|\__|\___/ ____/ -// Multi-Physics +// Multi-Physics // -// License: BSD License -// Kratos default license: kratos/license.txt +// License: BSD License +// Kratos default license: kratos/license.txt // // Main authors: Jordi Cotela Dalmau -// // -#if !defined(KRATOS_CARTESIAN_MESH_GENERATOR_H_INCLUDED ) -#define KRATOS_CARTESIAN_MESH_GENERATOR_H_INCLUDED +#pragma once // System includes -#include -#include // External includes // Project includes -#include "includes/define.h" +#include "containers/model.h" +#include "geometries/geometry.h" +#include "containers/pointer_vector.h" #include "modeler/modeler.h" -#include "spatial_containers/spatial_containers.h" namespace Kratos { -///@name Kratos Globals +///@addtogroup KratosCore ///@{ -///@} -///@name Type Definitions -///@{ - -///@} -///@name Enum's -///@{ - -///@} -///@name Functions -///@{ - -///@} ///@name Kratos Classes ///@{ -/// Short class definition. -/** Detail class definition. -*/ -class CartesianMeshGeneratorModeler : public Modeler +/** + * @class CartesianMeshGeneratorModeler + * @ingroup KratosCore + * @brief Generates a structured Cartesian mesh filling the bounding box of a + * source model part. + * @details For 3D problems each axis-aligned hexahedral cell is split into + * 6 tetrahedra using the Freudenthal (Kuhn) decomposition so that + * the resulting tetrahedral mesh is conforming across cell boundaries. + * For 2D problems quadrilateral elements are generated using a + * ray-casting inside test derived from the source boundary elements. + * + * The class is constructed with a reference source model part + * (from which the bounding box and boundary are taken) and a uniform + * element size. The generated nodes and elements are placed in the + * destination model part passed to GenerateCartesianMesh(). + * + * @author Jordi Cotela Dalmau + */ +class KRATOS_API(KRATOS_CORE) CartesianMeshGeneratorModeler + : public Modeler { public: ///@name Type Definitions @@ -58,562 +58,252 @@ class CartesianMeshGeneratorModeler : public Modeler /// Pointer definition of CartesianMeshGeneratorModeler KRATOS_CLASS_POINTER_DEFINITION(CartesianMeshGeneratorModeler); - typedef Modeler BaseType; - - typedef Point PointType; - - typedef Node NodeType; + /// The base class type + using BaseType = Modeler; - typedef Geometry GeometryType; + /// The nodes vector type definition + using NodesVectorType = PointerVector; - typedef PointerVector NodesVectorType; - - typedef std::size_t SizeType; + /// The size type definition + using SizeType = std::size_t; ///@} ///@name Life Cycle ///@{ - /// constructor. - CartesianMeshGeneratorModeler(ModelPart& rSourceModelPart, double ElementSize) : - mrModelPart(rSourceModelPart), mElementSize(ElementSize) - { - } + /** + * @brief Default constructor (required for registry). + * @note Do not call GenerateCartesianMesh() on a default-constructed instance. + */ + CartesianMeshGeneratorModeler(); + + /** + * @brief Constructor for use via the modeler factory / registry. + * @details Reads the input model part name, output model part name, + * element name and element size from @p ModelerParameters. + * The default parameters are: + * @code + * { + * "input_model_part_name" : "", + * "output_model_part_name" : "", + * "element_name" : "Element3D4N", + * "element_size" : 1.0 + * } + * @endcode + * SetupModelPart() performs the actual mesh generation. + * @param rModel Model that owns both the input and output model parts. + * @param ModelerParameters Parameters block (see default values above). + */ + CartesianMeshGeneratorModeler(Model& rModel, Parameters ModelerParameters = Parameters()); + + /** + * @brief Legacy direct constructor. + * @param rSourceModelPart Model part whose geometry defines the bounding + * box (and, for 2D, the boundary) used during mesh generation. + * @param ElementSize Uniform edge length of each voxel cell. The actual + * number of cells per direction is computed as + * ceil(bounding_box_extent / ElementSize). + */ + CartesianMeshGeneratorModeler(ModelPart& rSourceModelPart, double ElementSize) + : Modeler() + , mpModel(nullptr) + , mpSourceModelPart(&rSourceModelPart) + , mElementSize(ElementSize) {} /// Destructor. - virtual ~CartesianMeshGeneratorModeler() {} - - - ///@} - ///@name Operators - ///@{ - + ~CartesianMeshGeneratorModeler() override = default; ///@} ///@name Operations ///@{ - void GenerateMesh(ModelPart& rThisModelPart, Element const& rReferenceElement) - { - const unsigned int dimension = rReferenceElement.GetGeometry().Dimension(); - - KRATOS_WATCH(dimension); - - Timer::Start("Generating Mesh"); - - CalculateBoundingBox(mrModelPart, mMinPoint, mMaxPoint); - CalculateDivisionNumbers(); - - unsigned int start_node_id = mrModelPart.NumberOfNodes() + 1; - unsigned int start_element_id = 1; - //unsigned int segment_number_1 = mSegmentsNumber[0] + 1; - //unsigned int segment_number_2 = mSegmentsNumber[1] + 1; - //unsigned int segment_number_3 = mSegmentsNumber[2] + 1; - - - - const unsigned int number_of_nodes = mDivisionsNumber[0] * mDivisionsNumber[1] * mDivisionsNumber[2]; - //const unsigned int number_of_nodes = segment_number_1 * segment_number_2 * segment_number_3; - - const unsigned int number_of_elements = mSegmentsNumber[0] * mSegmentsNumber[1] * mSegmentsNumber[2]; - - KRATOS_WATCH(number_of_nodes); - - KRATOS_WATCH(number_of_elements); - - CalculateBoundaryIntersections(mrModelPart); - - CalculateIsInside(mrModelPart); - - CalculateNormals(); - - Timer::Start("Generating Nodes"); - - double x0 = mMinPoint.X(); - double y0 = mMinPoint.Y(); - double z0 = mMinPoint.Z(); - - ModelPart::NodesContainerType::ContainerType& nodes_array = rThisModelPart.NodesArray(); - ModelPart::NodesContainerType::ContainerType temp_nodes_array(number_of_nodes); - - ModelPart::ElementsContainerType::ContainerType& elements_array = rThisModelPart.ElementsArray(); - ModelPart::ElementsContainerType::ContainerType temp_elements_array(number_of_elements); - - for(unsigned int k = 0 ; k < mDivisionsNumber[2] ; k++) - for(unsigned int j = 0 ; j < mDivisionsNumber[1] ; j++) - for(unsigned int i = 0 ; i < mDivisionsNumber[0] ; i++) - temp_nodes_array[NodeIndex(i,j,k)] = NodeType::Pointer(new NodeType(start_node_id++, x0 + i * mElementSize, y0 + j * mElementSize, z0 + k * mElementSize)); - - unsigned int number_of_active_nodes = 0; - - for(unsigned int i = 0 ; i < number_of_nodes ; i++) - if(mIsInside[i]) - number_of_active_nodes++; - - nodes_array.resize(number_of_active_nodes); - - unsigned int index = 0; - - for(unsigned int i = 0 ; i < number_of_nodes ; i++) - if(mIsInside[i]) - nodes_array[index++]=temp_nodes_array[i]; - - Timer::Stop("Generating Nodes"); - - Timer::Start("Generating Elements"); - - Element::NodesArrayType element_nodes(4); - - if(dimension == 2) - { - unsigned int number_of_active_elements = 0; - for(unsigned int j = 0 ; j < mSegmentsNumber[1] ; j++) - for(unsigned int i = 0 ; i < mSegmentsNumber[0] ; i++) - { - if(mIsInside[NodeIndex(i,j,0)] & mIsInside[NodeIndex(i+1,j,0)] & mIsInside[NodeIndex(i+1,j+1,0)] & mIsInside[NodeIndex(i,j+1,0)]) - number_of_active_elements++; - } - - elements_array.resize(number_of_active_elements); - - unsigned int counter = 0; - - for(unsigned int j = 0 ; j < mSegmentsNumber[1] ; j++) - for(unsigned int i = 0 ; i < mSegmentsNumber[0] ; i++) - { - if(mIsInside[NodeIndex(i,j,0)] & mIsInside[NodeIndex(i+1,j,0)] & mIsInside[NodeIndex(i+1,j+1,0)] & mIsInside[NodeIndex(i,j+1,0)]) - { - element_nodes(0) = temp_nodes_array[NodeIndex(i,j,0)]; - element_nodes(1) = temp_nodes_array[NodeIndex(i+1,j,0)]; - element_nodes(2) = temp_nodes_array[NodeIndex(i+1,j+1,0)]; - element_nodes(3) = temp_nodes_array[NodeIndex(i,j+1,0)]; - - elements_array[counter++] = rReferenceElement.Create(start_element_id++, element_nodes, rReferenceElement.pGetProperties()); - //elements_array[ElementIndex(i,j,0)] = rReferenceElement.Create(start_element_id++, element_nodes, rReferenceElement.pGetProperties()); - } - } - for(ModelPart::ElementIterator i_element = mrModelPart.ElementsBegin() ; i_element != mrModelPart.ElementsEnd() ; i_element++) - { - Element::GeometryType& r_geometry = i_element->GetGeometry(); - - PointType point1 = r_geometry[0] + mNormals[r_geometry[0].Id()-1]; - PointType point2 = r_geometry[1] + mNormals[r_geometry[1].Id()-1]; - - unsigned int index1 = FindNearestNodeIndex(point1,mNormals[r_geometry[0].Id()-1]); - unsigned int index2 = FindNearestNodeIndex(point2,mNormals[r_geometry[1].Id()-1]); - - element_nodes(0) = r_geometry(0); - element_nodes(1) = r_geometry(1); - element_nodes(2) = temp_nodes_array[index2]; - element_nodes(3) = temp_nodes_array[index1]; - - elements_array.push_back(rReferenceElement.Create(start_element_id++, element_nodes, rReferenceElement.pGetProperties())); - } - } - - // TODO: assigning nodes and elements to the modelpart - // TODO: adding source boundary mesh nods to modelpart - - Timer::Stop("Generating Elements"); - - - - Timer::Stop("Generating Mesh"); - } - - - void CalculateNormals() - { - array_1d zero = ZeroVector(3); - - if(mNormals.size() != mrModelPart.NumberOfNodes()) - mNormals.resize(mrModelPart.NumberOfNodes(), zero); - else - std::fill(mNormals.begin(), mNormals.end(), zero); - - double coefficient = mElementSize / 2.00; - - for(ModelPart::ElementIterator i_element = mrModelPart.ElementsBegin() ; i_element != mrModelPart.ElementsEnd() ; i_element++) - { - Element::GeometryType& r_geometry = i_element->GetGeometry(); - array_1d normal; - - normal[0] = r_geometry[1].Y() - r_geometry[0].Y(); - normal[1] = - (r_geometry[1].X() - r_geometry[0].X()); - - normal *= coefficient / r_geometry.Length(); - - - mNormals[r_geometry[0].Id()-1] += normal ; - mNormals[r_geometry[1].Id()-1] += normal ; - } - for(ModelPart::NodeIterator i_node = mrModelPart.NodesBegin() ; i_node != mrModelPart.NodesEnd() ; i_node++) - noalias(i_node->FastGetSolutionStepValue(NORMAL)) = mNormals[i_node->Id()-1] ; - } - - unsigned int FindNearestNodeIndex(PointType& rThisPoint, array_1d& rNormal) - { - double x = (rThisPoint.X() - mMinPoint.X()) / mElementSize; - double y = (rThisPoint.Y() - mMinPoint.Y()) / mElementSize; - - unsigned int i = static_cast(x); - unsigned int j = static_cast(y); - - if(mIsInside[NodeIndex(i,j,0)]) - return NodeIndex(i,j,0); - - - if(rNormal[0] >= 0.00) - { - if(rNormal[1] >= 0.00) - { - if(rNormal[0] > rNormal[1]) - i++; - else - j++; - } - else - { - if(rNormal[0] > -rNormal[1]) - i++; - else - j--; - } - } - else - { - if(rNormal[1] >= 0.00) - { - if(-rNormal[0] > rNormal[1]) - i--; - else - j++; - } - else - { - if(-rNormal[0] > -rNormal[1]) - i--; - else - j--; - } - } - - - return NodeIndex(i,j,0); - - } - - void CalculateIsInside(ModelPart& rThisModelPart) - { - const unsigned int number_of_nodes = mDivisionsNumber[0] * mDivisionsNumber[1] * mDivisionsNumber[2]; - - if(mIsInside.size() != number_of_nodes) - { - mIsInside.resize(number_of_nodes, 0); - } - else - std::fill(mIsInside.begin(),mIsInside.end(), 0); - - int size = mSegmentsNumber[1] + 1; - - for(int j = 0 ; j < size ; j++) - { -// bool is_inside = false; - std::vector& j_intersections = mIntersections[j]; - - for(std::vector::iterator j_x = j_intersections.begin() ; j_x != j_intersections.end() ; j_x++) - { - std::vector::iterator start = j_x++; - - unsigned int i_start = static_cast(*start / mElementSize); - unsigned int i_end = static_cast(*j_x / mElementSize); - - for(unsigned int i = i_start + 1 ; i < i_end ; i++) - mIsInside[NodeIndex(i,j,0)]=true; - } - } - } - - void CalculateBoundaryIntersections(ModelPart& rThisModelPart) - { - std::vector number_of_intersections(mSegmentsNumber[1] + 1, 0); // number of intersections per row for mSegmentsNumber[1] + 1 rows in x direction - if(mIntersections.size() != mSegmentsNumber[1] + 1) - mIntersections.resize(mSegmentsNumber[1] + 1); // intersection points x coordinate per row for mSegmentsNumber[1] + 1 rows - - for(ModelPart::ElementIterator i_element = rThisModelPart.ElementsBegin() ; i_element != rThisModelPart.ElementsEnd() ; i_element++) - { - Element::GeometryType& r_geometry = i_element->GetGeometry(); - - double x1; - double x2; - double y1; - double y2; - if(r_geometry[0].Y() < r_geometry[1].Y()) - { - x1 = (r_geometry[0].X() - mMinPoint.X()); - x2 = (r_geometry[1].X() - mMinPoint.X()); - y1 = (r_geometry[0].Y() - mMinPoint.Y()); - y2 = (r_geometry[1].Y() - mMinPoint.Y()); - } - else - { - x1 = (r_geometry[1].X() - mMinPoint.X()); - x2 = (r_geometry[0].X() - mMinPoint.X()); - y1 = (r_geometry[1].Y() - mMinPoint.Y()); - y2 = (r_geometry[0].Y() - mMinPoint.Y()); - } - - unsigned int i_start = static_cast(y1 / mElementSize); - unsigned int i_end = static_cast(y2 / mElementSize); - - if (i_start*mElementSize < y1) - i_start++; - - double m = 0.00; - if(y1 != y2) - m = (x2 - x1) / (y2 - y1); - - double delta_x = m*mElementSize; - double x = x1 + (i_start * mElementSize - y1) * m ; - for(unsigned int i = i_start ; i <= i_end ; i++) - { - number_of_intersections[i]++; - mIntersections[i].push_back(x); - x += delta_x; - } - } - - for(int i = 0 ; i < mIntersections.size() ; i++) - std::sort(mIntersections[i].begin(), mIntersections[i].end()); - } - - void CalculateBoundingBox(ModelPart& rThisModelPart, Point& rMinPoint, Point& rMaxPoint) + /** + * @brief Creates a new instance via the modeler factory. + * @param rModel Owning model. + * @param ModelParameters Configuration parameters. + * @return Shared pointer to the new modeler. + */ + Modeler::Pointer Create(Model& rModel, const Parameters ModelParameters) const override { - if(rThisModelPart.NumberOfElements() == 0) - { - rMinPoint = PointType(); - rMaxPoint = PointType(); - return; - } - - if(rThisModelPart.ElementsBegin()->GetGeometry().empty()) - { - rMinPoint = PointType(); - rMaxPoint = PointType(); - return; - } - - - rMinPoint = rThisModelPart.ElementsBegin()->GetGeometry()[0]; - rMaxPoint = rMinPoint; - - for(ModelPart::ElementIterator i_element = rThisModelPart.ElementsBegin() ; - i_element != rThisModelPart.ElementsEnd() ; i_element++) - for(Element::GeometryType::iterator i_point = i_element->GetGeometry().begin() ; i_point != i_element->GetGeometry().end() ; i_point++) - for(unsigned int i = 0 ; i < PointType::Dimension() ; i++) - { - if(rMinPoint[i] > (*i_point)[i]) - rMinPoint[i] = (*i_point)[i]; - - if(rMaxPoint[i] < (*i_point)[i]) - rMaxPoint[i] = (*i_point)[i]; - } + return Kratos::make_shared(rModel, ModelParameters); } - - void CalculateDivisionNumbers() - { - if(mElementSize == 0.00) - return; - - for(unsigned int i = 0 ; i < PointType::Dimension() ; i++) - { - double delta = mMaxPoint[i] - mMinPoint[i]; - int segments_number = static_cast(delta / mElementSize); - - if (((segments_number * mElementSize) < delta)) - segments_number++; - - mSegmentsNumber[i] = segments_number; - KRATOS_WATCH(mSegmentsNumber[i]); - - mDivisionsNumber[i] = segments_number + 1; - - if(mSegmentsNumber[i] == 0) - mSegmentsNumber[i]++; - } - - } - - virtual void GenerateNodes(ModelPart& ThisModelPart) - { - //std::vector - } - - + /** + * @brief Returns the default parameters accepted by the Model constructor. + * @return A Parameters object with default values. + */ + const Parameters GetDefaultParameters() const override; + + /** + * @brief Runs mesh generation when the modeler is used through the factory. + * @details Reads the input/output model part names and element name from + * the stored parameters and delegates to GenerateCartesianMesh(). + */ + void SetupModelPart() override; + + /** + * @brief Generates a Cartesian mesh in the destination model part. + * @details The element type (and therefore the dimension) is determined + * by @p rReferenceElement. For 3D elements six tetrahedra are + * created per hex cell (Freudenthal decomposition). For 2D + * elements one quadrilateral is created per active cell. + * @param rThisModelPart Destination model part that will receive the + * generated nodes and elements. + * @param rReferenceElement Prototype element whose geometry dimension + * selects the 2D or 3D code path and whose properties are + * reused for every created element. + */ + void GenerateCartesianMesh(ModelPart& rThisModelPart, Element const& rReferenceElement); + + /** + * @brief Computes outward-pointing normals at every boundary node of + * the source model part (2D only). + * @details The normals are scaled to half the element size and + * accumulated over all boundary elements sharing each node. + * The result is stored in the NORMAL solution-step variable + * and in the internal @c mNormals cache used by + * FindNearestNodeIndex(). + */ + void CalculateNormals(); + + /** + * @brief Returns the flat index of the grid node closest to @p rThisPoint + * in the direction given by @p rNormal (2D only). + * @details Starting from the cell that contains the point the method + * steps one cell in the dominant normal direction until it + * finds a node that is marked as inside. + * @param rThisPoint Query point in world coordinates. + * @param rNormal Outward boundary normal at the query point; used to + * determine which neighbouring cell to prefer. + * @return Flat index into the internal node array of the nearest inside + * node. + */ + unsigned int FindNearestNodeIndex(Point& rThisPoint, array_1d& rNormal); + + /** + * @brief Marks each grid node as inside or outside the 2D boundary. + * @details Uses the precomputed horizontal ray–edge intersection table + * (populated by CalculateBoundaryIntersections()) to flood-fill + * the grid cells between pairs of intersections on each + * horizontal scan line. + * @param rThisModelPart The source model part providing the boundary. + */ + void CalculateIsInside(ModelPart& rThisModelPart); + + /** + * @brief Builds the horizontal ray–edge intersection table (2D only). + * @details For every boundary element in @p rThisModelPart the method + * computes the x-coordinates at which the element edge crosses + * each horizontal grid line and stores them in @c mIntersections. + * Each row of the table is sorted in ascending order after all + * edges have been processed. + * @param rThisModelPart The source model part providing the boundary edges. + */ + void CalculateBoundaryIntersections(ModelPart& rThisModelPart); + + /** + * @brief Computes the axis-aligned bounding box of @p rThisModelPart. + * @details When the model part contains elements the bounding box is + * computed from element vertices. When it contains only nodes + * (e.g. a condition-only mesh read from an STL file) the bounding + * box is computed directly from the node coordinates. An empty + * model part yields the origin for both corners. + * @param rThisModelPart The model part to bound. + * @param rMinPoint Output: minimum corner of the bounding box. + * @param rMaxPoint Output: maximum corner of the bounding box. + */ + void CalculateBoundingBox(ModelPart& rThisModelPart, Point& rMinPoint, Point& rMaxPoint); + + /** + * @brief Computes the number of grid cells and node layers per direction. + * @details Sets @c mSegmentsNumber[i] = ceil(extent_i / mElementSize) + * and @c mDivisionsNumber[i] = mSegmentsNumber[i] + 1. + * Directions with zero extent are given one segment so that the + * degenerate case does not crash. + */ + void CalculateDivisionNumbers(); + + /** + * @brief Returns the flat element index for cell (i, j, k). + * @param i Cell index along X. + * @param j Cell index along Y. + * @param k Cell index along Z. + * @return Flat index = i + nx*j + nx*ny*k, where nx = mDivisionsNumber[0] + * and ny = mDivisionsNumber[1]. + */ unsigned int ElementIndex(unsigned int i, unsigned int j, unsigned int k) { return i + mDivisionsNumber[0] * j + mDivisionsNumber[0] * mDivisionsNumber[1] * k; } + /** + * @brief Returns the flat node index for grid point (i, j, k). + * @param i Node index along X. + * @param j Node index along Y. + * @param k Node index along Z. + * @return Flat index = i + sx*j + sx*sy*k, where sx = mSegmentsNumber[0] + * and sy = mSegmentsNumber[1]. + */ unsigned int NodeIndex(unsigned int i, unsigned int j, unsigned int k) { return i + mSegmentsNumber[0] * j + mSegmentsNumber[0] * mSegmentsNumber[1] * k; } - - - - - ///@} - ///@name Access - ///@{ - - - ///@} - ///@name Inquiry - ///@{ - - ///@} ///@name Input and output ///@{ - /// Turn back information as a string. - virtual std::string Info() const - { - return "CartesianMeshGeneratorModeler"; - } - - /// Print information about this object. - virtual void PrintInfo(std::ostream& rOStream) const - { - rOStream << Info(); - } - - /// Print object's data. - virtual void PrintData(std::ostream& rOStream) const - { - } - - - ///@} - ///@name Friends - ///@{ - - - ///@} - -protected: - ///@name Protected static Member Variables - ///@{ - - - ///@} - ///@name Protected member Variables - ///@{ - + /// @brief Returns the class name as a string. + std::string Info() const override; - ///@} - ///@name Protected Operators - ///@{ - - - ///@} - ///@name Protected Operations - ///@{ - - - ///@} - ///@name Protected Access - ///@{ - - - ///@} - ///@name Protected Inquiry - ///@{ - - - ///@} - ///@name Protected LifeCycle - ///@{ + /// @brief Prints the class name to @p rOStream. + void PrintInfo(std::ostream& rOStream) const override; + /// @brief Prints object data to @p rOStream (currently empty). + void PrintData(std::ostream& rOStream) const override; ///@} private: - ///@name Static Member Variables - ///@{ - - - ///@} ///@name Member Variables ///@{ - ModelPart& mrModelPart; - - double mElementSize; - PointType mMinPoint; - PointType mMaxPoint; - - unsigned int mSegmentsNumber[3]; + /// Owning model (only set when using the Model constructor). + Model* mpModel = nullptr; - unsigned int mDivisionsNumber[3]; + /// Source model part whose geometry defines the domain to mesh (nullptr only for the default-constructed prototype). + ModelPart* mpSourceModelPart; - std::vector > mIntersections; - std::vector mIsInside; + /// Uniform voxel cell edge length. + double mElementSize = 0.0; - std::vector > mNormals; + /// Minimum corner of the bounding box of the source model part. + Point mMinPoint; - ///@} - ///@name Private Operators - ///@{ + /// Maximum corner of the bounding box of the source model part. + Point mMaxPoint; + /// Number of voxel cells per direction (X=0, Y=1, Z=2). + unsigned int mSegmentsNumber[3] = {0, 0, 0}; - ///@} - ///@name Private Operations - ///@{ + /// Number of node layers per direction: mSegmentsNumber[i] + 1. + unsigned int mDivisionsNumber[3] = {0, 0, 0}; + /// Sorted x-coordinates of edge–scanline intersections per row (2D only). + std::vector> mIntersections; - void GenerateNodes(ModelPart& ThisModelPart, GeometryType& rGeometry, SizeType NumberOfSegments, SizeType StartNodeId) - { - double x1 = rGeometry[0][0]; - double y1 = rGeometry[0][1]; - double z1 = rGeometry[0][2]; - double x2 = rGeometry[1][0]; - double y2 = rGeometry[1][1]; - double z2 = rGeometry[1][2]; - - double dx = (x2 - x1) / NumberOfSegments; - double dy = (y2 - y1) / NumberOfSegments; - double dz = (z2 - z1) / NumberOfSegments; - - for(SizeType i = 1 ; i < NumberOfSegments - 1 ; i++) - { - ThisModelPart.CreateNewNode(StartNodeId++, x1 + i * dx, y1 + i * dy, z1 + i * dz); - } - } - - ///@} - ///@name Private Access - ///@{ - - - ///@} - ///@name Private Inquiry - ///@{ + /// Per-node inside flag populated by CalculateIsInside() (2D only). + std::vector mIsInside; + /// Accumulated boundary normals at each source-model-part node (2D only). + std::vector> mNormals; ///@} - ///@name Un accessible methods + ///@name Inaccessible methods ///@{ - /// Assignment operator. - CartesianMeshGeneratorModeler& operator=(CartesianMeshGeneratorModeler const& rOther); - - /// Copy constructor. - CartesianMeshGeneratorModeler(CartesianMeshGeneratorModeler const& rOther); - + CartesianMeshGeneratorModeler& operator=(CartesianMeshGeneratorModeler const& rOther) = delete; + CartesianMeshGeneratorModeler(CartesianMeshGeneratorModeler const& rOther) = delete; ///@} @@ -621,34 +311,17 @@ class CartesianMeshGeneratorModeler : public Modeler ///@} -///@name Type Definitions -///@{ +/// @brief Input stream operator (no-op). +inline std::istream& operator>>(std::istream& rIStream, CartesianMeshGeneratorModeler& rThis) { return rIStream; } - -///@} -///@name Input and output -///@{ - - -/// input stream function -inline std::istream& operator >> (std::istream& rIStream, - CartesianMeshGeneratorModeler& rThis); - -/// output stream function -inline std::ostream& operator << (std::ostream& rOStream, - const CartesianMeshGeneratorModeler& rThis) +/// @brief Output stream operator — prints Info() followed by PrintData(). +inline std::ostream& operator<<(std::ostream& rOStream, const CartesianMeshGeneratorModeler& rThis) { rThis.PrintInfo(rOStream); rOStream << std::endl; rThis.PrintData(rOStream); - return rOStream; } -///@} - - -} // namespace Kratos. - -#endif // KRATOS_CARTESIAN_MESH_GENERATOR_H_INCLUDED defined +} // namespace Kratos diff --git a/kratos/python/add_modeler_to_python.cpp b/kratos/python/add_modeler_to_python.cpp index 45e0f56f86b2..d39a2c01deda 100644 --- a/kratos/python/add_modeler_to_python.cpp +++ b/kratos/python/add_modeler_to_python.cpp @@ -28,6 +28,7 @@ #include "modeler/combine_model_part_modeler.h" #include "modeler/voxel_mesh_generator_modeler.h" #include "modeler/clean_up_problematic_triangles_modeler.h" +#include "modeler/cartesian_mesh_generator_modeler.h" namespace Kratos::Python { @@ -122,6 +123,14 @@ void AddModelerToPython(pybind11::module& m) py::class_(m, "CleanUpProblematicTrianglesModeler") .def(py::init()) ; + + py::class_(m, "CartesianMeshGeneratorModeler") + .def(py::init()) + .def(py::init()) + .def("GenerateMesh", [](CartesianMeshGeneratorModeler& self, ModelPart& rModelPart, const std::string& rElementName) { + self.GenerateCartesianMesh(rModelPart, KratosComponents::Get(rElementName)); + }) + ; } } // namespace Kratos::Python. \ No newline at end of file diff --git a/kratos/sources/kratos_application.cpp b/kratos/sources/kratos_application.cpp index 1cd550418682..afa802ea5b1a 100644 --- a/kratos/sources/kratos_application.cpp +++ b/kratos/sources/kratos_application.cpp @@ -272,6 +272,7 @@ void KratosApplication::RegisterKratosCore() { KRATOS_REGISTER_MODELER("CombineModelPartModeler", mCombineModelPartModeler); KRATOS_REGISTER_MODELER("ConnectivityPreserveModeler", mConnectivityPreserveModeler); KRATOS_REGISTER_MODELER("VoxelMeshGeneratorModeler", mVoxelMeshGeneratorModeler); + KRATOS_REGISTER_MODELER("CartesianMeshGeneratorModeler", mCartesianMeshGeneratorModeler); KRATOS_REGISTER_MODELER("CleanUpProblematicTrianglesModeler", mCleanUpProblematicTrianglesModeler); // Register general geometries: diff --git a/kratos/tests/cpp_tests/modeler/test_cartesian_mesh_generator_modeler.cpp b/kratos/tests/cpp_tests/modeler/test_cartesian_mesh_generator_modeler.cpp new file mode 100644 index 000000000000..5823a0304d14 --- /dev/null +++ b/kratos/tests/cpp_tests/modeler/test_cartesian_mesh_generator_modeler.cpp @@ -0,0 +1,192 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Jordi Cotela Dalmau +// + +// System includes + +// External includes + +// Project includes +#include "testing/testing.h" +#include "containers/model.h" +#include "includes/element.h" +#include "includes/kratos_components.h" +#include "modeler/cartesian_mesh_generator_modeler.h" + +namespace Kratos::Testing +{ + +namespace { + +/// Build a square 2x2 bounding-box model part with 8 corner nodes only (no elements), +/// so the node-based BB fallback is exercised. +ModelPart& CreateBoxSourceModelPart(Model& rModel, const std::string& rName, + double xMin, double yMin, double zMin, + double xMax, double yMax, double zMax) +{ + auto& r_mp = rModel.CreateModelPart(rName); + r_mp.CreateNewNode(1, xMin, yMin, zMin); + r_mp.CreateNewNode(2, xMax, yMin, zMin); + r_mp.CreateNewNode(3, xMax, yMax, zMin); + r_mp.CreateNewNode(4, xMin, yMax, zMin); + r_mp.CreateNewNode(5, xMin, yMin, zMax); + r_mp.CreateNewNode(6, xMax, yMin, zMax); + r_mp.CreateNewNode(7, xMax, yMax, zMax); + r_mp.CreateNewNode(8, xMin, yMax, zMax); + return r_mp; +} + +} // anonymous namespace + +/***********************************************************************************/ +/***********************************************************************************/ + +/// Test that GenerateCartesianMesh (direct API) fills a 3D model part with the correct +/// node and element counts for a unit-cube source with a given element size. +KRATOS_TEST_CASE_IN_SUITE(CartesianMeshGeneratorModeler3DNodeCount, KratosCoreFastSuite) +{ + Model model; + // 1x1x1 box, element size 0.5 → 2 segments per direction → 3^3=27 nodes, 2^3*6=48 tets + auto& r_source = CreateBoxSourceModelPart(model, "source", 0.0, 0.0, 0.0, 1.0, 1.0, 1.0); + auto& r_output = model.CreateModelPart("output"); + r_output.CreateNewProperties(0); + + CartesianMeshGeneratorModeler modeler(r_source, 0.5); + modeler.GenerateCartesianMesh(r_output, KratosComponents::Get("Element3D4N")); + + // 3 nodes per direction → 27 total + KRATOS_EXPECT_EQ(r_output.NumberOfNodes(), 27u); + // 2 cells per direction → 8 hex cells × 6 tets = 48 tetrahedra + KRATOS_EXPECT_EQ(r_output.NumberOfElements(), 48u); +} + +/***********************************************************************************/ +/***********************************************************************************/ + +/// Test that every generated node lies inside (or on the boundary of) the source BB. +KRATOS_TEST_CASE_IN_SUITE(CartesianMeshGeneratorModeler3DNodesInsideBB, KratosCoreFastSuite) +{ + Model model; + auto& r_source = CreateBoxSourceModelPart(model, "source", -1.0, -1.0, -1.0, 1.0, 1.0, 1.0); + auto& r_output = model.CreateModelPart("output"); + r_output.CreateNewProperties(0); + + CartesianMeshGeneratorModeler modeler(r_source, 0.5); + modeler.GenerateCartesianMesh(r_output, KratosComponents::Get("Element3D4N")); + + for (const auto& r_node : r_output.Nodes()) { + KRATOS_EXPECT_GE(r_node.X(), -1.0 - 1e-12); + KRATOS_EXPECT_LE(r_node.X(), 1.0 + 1e-12); + KRATOS_EXPECT_GE(r_node.Y(), -1.0 - 1e-12); + KRATOS_EXPECT_LE(r_node.Y(), 1.0 + 1e-12); + KRATOS_EXPECT_GE(r_node.Z(), -1.0 - 1e-12); + KRATOS_EXPECT_LE(r_node.Z(), 1.0 + 1e-12); + } +} + +/***********************************************************************************/ +/***********************************************************************************/ + +/// Test that each element has exactly 4 nodes (tetrahedron) and that all node IDs +/// reference nodes that are present in the output model part. +KRATOS_TEST_CASE_IN_SUITE(CartesianMeshGeneratorModeler3DElementConnectivity, KratosCoreFastSuite) +{ + Model model; + auto& r_source = CreateBoxSourceModelPart(model, "source", 0.0, 0.0, 0.0, 1.0, 1.0, 1.0); + auto& r_output = model.CreateModelPart("output"); + r_output.CreateNewProperties(0); + + CartesianMeshGeneratorModeler modeler(r_source, 1.0); + modeler.GenerateCartesianMesh(r_output, KratosComponents::Get("Element3D4N")); + + // 1 cell per direction, 6 tets + KRATOS_EXPECT_EQ(r_output.NumberOfElements(), 6u); + + for (const auto& r_elem : r_output.Elements()) { + KRATOS_EXPECT_EQ(r_elem.GetGeometry().size(), 4u); + for (unsigned int n = 0; n < 4; ++n) { + const auto node_id = r_elem.GetGeometry()[n].Id(); + KRATOS_EXPECT_TRUE(r_output.HasNode(node_id)); + } + } +} + +/***********************************************************************************/ +/***********************************************************************************/ + +/// Test the factory / registry constructor (Model + Parameters) via SetupModelPart(). +KRATOS_TEST_CASE_IN_SUITE(CartesianMeshGeneratorModelerSetupModelPart, KratosCoreFastSuite) +{ + Model model; + CreateBoxSourceModelPart(model, "source", 0.0, 0.0, 0.0, 2.0, 2.0, 2.0); + + Parameters params(R"({ + "input_model_part_name" : "source", + "output_model_part_name" : "output", + "element_name" : "Element3D4N", + "element_size" : 1.0 + })"); + + CartesianMeshGeneratorModeler modeler(model, params); + modeler.SetupModelPart(); + + const auto& r_output = model.GetModelPart("output"); + // 2x2x2 box, element size 1 → 3^3=27 nodes, 2^3*6=48 tets + KRATOS_EXPECT_EQ(r_output.NumberOfNodes(), 27u); + KRATOS_EXPECT_EQ(r_output.NumberOfElements(), 48u); +} + +/***********************************************************************************/ +/***********************************************************************************/ + +/// Test that the bounding box helper correctly handles a node-only model part +/// (no elements), which is the typical output of StlIO (conditions only). +KRATOS_TEST_CASE_IN_SUITE(CartesianMeshGeneratorModelerBoundingBoxNodeOnly, KratosCoreFastSuite) +{ + Model model; + auto& r_source = CreateBoxSourceModelPart(model, "source", -2.0, -3.0, -4.0, 2.0, 3.0, 4.0); + auto& r_output = model.CreateModelPart("output"); + r_output.CreateNewProperties(0); + + CartesianMeshGeneratorModeler modeler(r_source, 2.0); + Point min_pt, max_pt; + modeler.CalculateBoundingBox(r_source, min_pt, max_pt); + + KRATOS_EXPECT_NEAR(min_pt.X(), -2.0, 1e-12); + KRATOS_EXPECT_NEAR(min_pt.Y(), -3.0, 1e-12); + KRATOS_EXPECT_NEAR(min_pt.Z(), -4.0, 1e-12); + KRATOS_EXPECT_NEAR(max_pt.X(), 2.0, 1e-12); + KRATOS_EXPECT_NEAR(max_pt.Y(), 3.0, 1e-12); + KRATOS_EXPECT_NEAR(max_pt.Z(), 4.0, 1e-12); +} + +/***********************************************************************************/ +/***********************************************************************************/ + +/// Test CalculateDivisionNumbers: exact multiples produce the right segment count. +KRATOS_TEST_CASE_IN_SUITE(CartesianMeshGeneratorModelerDivisionNumbers, KratosCoreFastSuite) +{ + Model model; + // 3x2x1 box, element size 0.5 → 6, 4, 2 segments + auto& r_source = CreateBoxSourceModelPart(model, "source", 0.0, 0.0, 0.0, 3.0, 2.0, 1.0); + auto& r_output = model.CreateModelPart("output"); + r_output.CreateNewProperties(0); + + CartesianMeshGeneratorModeler modeler(r_source, 0.5); + modeler.GenerateCartesianMesh(r_output, KratosComponents::Get("Element3D4N")); + + // (6+1)*(4+1)*(2+1) = 7*5*3 = 105 nodes + KRATOS_EXPECT_EQ(r_output.NumberOfNodes(), 105u); + // 6*4*2*6 = 288 tetrahedra + KRATOS_EXPECT_EQ(r_output.NumberOfElements(), 288u); +} + +} // namespace Kratos::Testing diff --git a/kratos/tests/test_KratosCore.py b/kratos/tests/test_KratosCore.py index 94057103c98e..78333d9cc85a 100644 --- a/kratos/tests/test_KratosCore.py +++ b/kratos/tests/test_KratosCore.py @@ -100,6 +100,7 @@ import test_import_obj_modeler import test_vectorized_interpolation import test_clean_up_problematic_triangles_modeler +import test_cartesian_mesh_generator_modeler import test_tetrahedral_mesh_orientation_check import test_nd_data import test_tensor_adaptors @@ -108,6 +109,7 @@ import test_constraint_restart import test_vtu_output import test_ensight_output_process +import test_vtu_output # Import modules required for sequential orchestrator test from test_sequential_orchestrator import EmptyAnalysisStage @@ -231,6 +233,7 @@ def AssembleTestSuites(): smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_import_obj_modeler.TestImportOBJModeler])) smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_vectorized_interpolation.TestVectorizedInterpolation])) smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_clean_up_problematic_triangles_modeler.TestCleanUpProblematicTrianglesModeler])) + smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_cartesian_mesh_generator_modeler.TestCartesianMeshGeneratorModeler])) smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_tetrahedral_mesh_orientation_check.TestTetrahedralMeshOrientationCheck])) smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_nd_data.TestNDData])) smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_tensor_adaptors.TestTensorAdaptors])) @@ -241,6 +244,9 @@ def AssembleTestSuites(): smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_vtu_output.TestVtuOutput2D])) smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_vtu_output.TestVtuOutput3D])) smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_ensight_output_process.TestEnsightOutputProcess])) + smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_vtu_output.TestVtuOutput])) + smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_vtu_output.TestVtuOutput2D])) + smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_vtu_output.TestVtuOutput3D])) if sympy_available: smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_sympy_fe_utilities.TestSympyFEUtilities])) diff --git a/kratos/tests/test_cartesian_mesh_generator_modeler.py b/kratos/tests/test_cartesian_mesh_generator_modeler.py new file mode 100644 index 000000000000..6b4495aeed86 --- /dev/null +++ b/kratos/tests/test_cartesian_mesh_generator_modeler.py @@ -0,0 +1,188 @@ +import KratosMultiphysics +import KratosMultiphysics.KratosUnittest as KratosUnittest + + +def _CreateBoxSourceModelPart(model, name, x_min, y_min, z_min, x_max, y_max, z_max): + """Create a model part containing only the 8 corner nodes of an axis-aligned box.""" + mp = model.CreateModelPart(name) + mp.CreateNewNode(1, x_min, y_min, z_min) + mp.CreateNewNode(2, x_max, y_min, z_min) + mp.CreateNewNode(3, x_max, y_max, z_min) + mp.CreateNewNode(4, x_min, y_max, z_min) + mp.CreateNewNode(5, x_min, y_min, z_max) + mp.CreateNewNode(6, x_max, y_min, z_max) + mp.CreateNewNode(7, x_max, y_max, z_max) + mp.CreateNewNode(8, x_min, y_max, z_max) + return mp + + +class TestCartesianMeshGeneratorModeler(KratosUnittest.TestCase): + + # ------------------------------------------------------------------ + # Direct API (ModelPart + element_size constructor) + # ------------------------------------------------------------------ + + def test_3d_node_count(self): + """1x1x1 box, element_size=0.5 → 3^3=27 nodes, 2^3*6=48 tets.""" + model = KratosMultiphysics.Model() + source = _CreateBoxSourceModelPart(model, "source", 0.0, 0.0, 0.0, 1.0, 1.0, 1.0) + output = model.CreateModelPart("output") + output.CreateNewProperties(0) + + modeler = KratosMultiphysics.CartesianMeshGeneratorModeler(source, 0.5) + modeler.GenerateMesh(output, "Element3D4N") + + self.assertEqual(output.NumberOfNodes(), 27) + self.assertEqual(output.NumberOfElements(), 48) + + def test_3d_nodes_inside_bounding_box(self): + """All generated nodes must lie inside (or on the boundary of) the source BB.""" + model = KratosMultiphysics.Model() + source = _CreateBoxSourceModelPart(model, "source", -1.0, -1.0, -1.0, 1.0, 1.0, 1.0) + output = model.CreateModelPart("output") + output.CreateNewProperties(0) + + modeler = KratosMultiphysics.CartesianMeshGeneratorModeler(source, 0.5) + modeler.GenerateMesh(output, "Element3D4N") + + tol = 1e-12 + for node in output.Nodes: + self.assertGreaterEqual(node.X, -1.0 - tol) + self.assertLessEqual(node.X, 1.0 + tol) + self.assertGreaterEqual(node.Y, -1.0 - tol) + self.assertLessEqual(node.Y, 1.0 + tol) + self.assertGreaterEqual(node.Z, -1.0 - tol) + self.assertLessEqual(node.Z, 1.0 + tol) + + def test_3d_element_connectivity(self): + """Each element must have 4 nodes that exist in the output model part.""" + model = KratosMultiphysics.Model() + source = _CreateBoxSourceModelPart(model, "source", 0.0, 0.0, 0.0, 1.0, 1.0, 1.0) + output = model.CreateModelPart("output") + output.CreateNewProperties(0) + + # element_size = 1.0 → 1 cell per direction → 6 tetrahedra + modeler = KratosMultiphysics.CartesianMeshGeneratorModeler(source, 1.0) + modeler.GenerateMesh(output, "Element3D4N") + + self.assertEqual(output.NumberOfElements(), 6) + node_ids = {node.Id for node in output.Nodes} + for elem in output.Elements: + self.assertEqual(len(elem.GetNodes()), 4) + for node in elem.GetNodes(): + self.assertIn(node.Id, node_ids) + + def test_3d_non_unit_box(self): + """3x2x1 box, element_size=0.5 → 7*5*3=105 nodes, 6*4*2*6=288 tets.""" + model = KratosMultiphysics.Model() + source = _CreateBoxSourceModelPart(model, "source", 0.0, 0.0, 0.0, 3.0, 2.0, 1.0) + output = model.CreateModelPart("output") + output.CreateNewProperties(0) + + modeler = KratosMultiphysics.CartesianMeshGeneratorModeler(source, 0.5) + modeler.GenerateMesh(output, "Element3D4N") + + self.assertEqual(output.NumberOfNodes(), 105) + self.assertEqual(output.NumberOfElements(), 288) + + def test_3d_uniform_spacing(self): + """Node coordinates should be uniformly spaced by element_size.""" + model = KratosMultiphysics.Model() + source = _CreateBoxSourceModelPart(model, "source", 0.0, 0.0, 0.0, 1.0, 1.0, 1.0) + output = model.CreateModelPart("output") + output.CreateNewProperties(0) + + element_size = 0.5 + modeler = KratosMultiphysics.CartesianMeshGeneratorModeler(source, element_size) + modeler.GenerateMesh(output, "Element3D4N") + + # Collect unique X values and check spacing + xs = sorted({round(node.X, 10) for node in output.Nodes}) + for i in range(1, len(xs)): + self.assertAlmostEqual(xs[i] - xs[i-1], element_size, places=10) + + # ------------------------------------------------------------------ + # Factory / registry constructor (Model + Parameters) + # ------------------------------------------------------------------ + + def test_setup_model_part(self): + """SetupModelPart() via factory constructor must produce the same mesh.""" + model = KratosMultiphysics.Model() + _CreateBoxSourceModelPart(model, "source", 0.0, 0.0, 0.0, 2.0, 2.0, 2.0) + + params = KratosMultiphysics.Parameters("""{ + "input_model_part_name" : "source", + "output_model_part_name" : "output", + "element_name" : "Element3D4N", + "element_size" : 1.0 + }""") + + modeler = KratosMultiphysics.CartesianMeshGeneratorModeler(model, params) + modeler.SetupModelPart() + + output = model.GetModelPart("output") + # 2x2x2 box, element_size=1 → 3^3=27 nodes, 2^3*6=48 tets + self.assertEqual(output.NumberOfNodes(), 27) + self.assertEqual(output.NumberOfElements(), 48) + + def test_setup_model_part_creates_output(self): + """SetupModelPart() must create the output model part if it does not exist.""" + model = KratosMultiphysics.Model() + _CreateBoxSourceModelPart(model, "box", 0.0, 0.0, 0.0, 1.0, 1.0, 1.0) + + params = KratosMultiphysics.Parameters("""{ + "input_model_part_name" : "box", + "output_model_part_name" : "mesh", + "element_name" : "Element3D4N", + "element_size" : 1.0 + }""") + + KratosMultiphysics.CartesianMeshGeneratorModeler(model, params).SetupModelPart() + self.assertTrue(model.HasModelPart("mesh")) + + def test_missing_input_name_raises(self): + """Empty input_model_part_name must raise an error.""" + model = KratosMultiphysics.Model() + model.CreateModelPart("source") + + params = KratosMultiphysics.Parameters("""{ + "input_model_part_name" : "", + "output_model_part_name" : "output", + "element_name" : "Element3D4N", + "element_size" : 1.0 + }""") + + modeler = KratosMultiphysics.CartesianMeshGeneratorModeler(model, params) + with self.assertRaises(RuntimeError): + modeler.SetupModelPart() + + def test_missing_output_name_raises(self): + """Empty output_model_part_name must raise an error.""" + model = KratosMultiphysics.Model() + _CreateBoxSourceModelPart(model, "source", 0.0, 0.0, 0.0, 1.0, 1.0, 1.0) + + params = KratosMultiphysics.Parameters("""{ + "input_model_part_name" : "source", + "output_model_part_name" : "", + "element_name" : "Element3D4N", + "element_size" : 1.0 + }""") + + modeler = KratosMultiphysics.CartesianMeshGeneratorModeler(model, params) + with self.assertRaises(RuntimeError): + modeler.SetupModelPart() + + # ------------------------------------------------------------------ + # Info / string interface + # ------------------------------------------------------------------ + + def test_info(self): + """Info() must return the class name.""" + model = KratosMultiphysics.Model() + source = model.CreateModelPart("source") + modeler = KratosMultiphysics.CartesianMeshGeneratorModeler(source, 1.0) + self.assertEqual(str(modeler).strip(), "CartesianMeshGeneratorModeler") + +if __name__ == '__main__': + KratosMultiphysics.Logger.GetDefaultOutput().SetSeverity(KratosMultiphysics.Logger.Severity.WARNING) + KratosUnittest.main()