Skip to content

Commit

Permalink
Merge remote-tracking branch 'origin/master' into optapp/windows_comp…
Browse files Browse the repository at this point in the history 
…ilation
  • Loading branch information
@sunethwarna
sunethwarna committed Aug 3, 2023
2 parents 495b3a9 + de6dff8 commit ca793fd
Show file tree
Hide file tree
Showing 99 changed files with 5,646 additions and 1,202 deletions.
131 changes: 78 additions & 53 deletions applications/IgaApplication/custom_modelers/nurbs_geometry_modeler.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,22 +22,39 @@ namespace Kratos
///@{

void NurbsGeometryModeler::SetupGeometryModel(){
// Get bounding box
KRATOS_ERROR_IF_NOT(mParameters.Has("lower_point"))
<< "NurbsGeometryModeler: Missing \"lower_point\" section" << std::endl;
// Get bounding box physical space.
KRATOS_ERROR_IF_NOT(mParameters.Has("lower_point_xyz"))
<< "NurbsGeometryModeler: Missing \"lower_point_xyz\" section" << std::endl;

KRATOS_ERROR_IF_NOT( mParameters["lower_point"].GetVector().size() == 3 )
<< "NurbsGeometryModeler: \"lower_point\" must be defined in 3 dimensions." << std::endl;
KRATOS_ERROR_IF_NOT( mParameters["lower_point_xyz"].GetVector().size() == 3 )
<< "NurbsGeometryModeler: \"lower_point_xyz\" must be defined in 3 dimensions." << std::endl;

Point point_a(mParameters["lower_point"].GetVector());
Point point_a_xyz(mParameters["lower_point_xyz"].GetVector());

KRATOS_ERROR_IF_NOT(mParameters.Has("upper_point"))
<< "NurbsGeometryModeler: Missing \"upper_point\" section" << std::endl;
KRATOS_ERROR_IF_NOT(mParameters.Has("upper_point_xyz"))
<< "NurbsGeometryModeler: Missing \"upper_point_xyz\" section" << std::endl;

KRATOS_ERROR_IF_NOT( mParameters["upper_point"].GetVector().size() == 3 )
<< "NurbsGeometryModeler: \"upper_point\" must be defined in 3 dimensions." << std::endl;
KRATOS_ERROR_IF_NOT( mParameters["upper_point_xyz"].GetVector().size() == 3 )
<< "NurbsGeometryModeler: \"upper_point_xyz\" must be defined in 3 dimensions." << std::endl;

Point point_b(mParameters["upper_point"].GetVector());
Point point_b_xyz(mParameters["upper_point_xyz"].GetVector());

// Get bounding box parametric space.
KRATOS_ERROR_IF_NOT(mParameters.Has("lower_point_uvw"))
<< "NurbsGeometryModeler: Missing \"lower_point_uvw\" section" << std::endl;

KRATOS_ERROR_IF_NOT( mParameters["lower_point_uvw"].GetVector().size() == 3 )
<< "NurbsGeometryModeler: \"lower_point_uvw\" must be defined in 3 dimensions." << std::endl;

Point point_a_uvw(mParameters["lower_point_uvw"].GetVector());

KRATOS_ERROR_IF_NOT(mParameters.Has("upper_point_uvw"))
<< "NurbsGeometryModeler: Missing \"upper_point_uvw\" section" << std::endl;

KRATOS_ERROR_IF_NOT( mParameters["upper_point_uvw"].GetVector().size() == 3 )
<< "NurbsGeometryModeler: \"upper_point_uvw\" must be defined in 3 dimensions." << std::endl;

Point point_b_uvw(mParameters["upper_point_uvw"].GetVector());

// Get polynomial order and number of knotspans
KRATOS_ERROR_IF_NOT(mParameters.Has("polynomial_order"))
Expand Down Expand Up @@ -67,7 +84,7 @@ namespace Kratos
SizeType num_knot_span_u = mParameters["number_of_knot_spans"].GetArrayItem(0).GetInt();
SizeType num_knot_span_v = mParameters["number_of_knot_spans"].GetArrayItem(1).GetInt();

CreateAndAddRegularGrid2D(model_part, point_a, point_b, p_u, p_v, num_knot_span_u, num_knot_span_v);
CreateAndAddRegularGrid2D(model_part, point_a_xyz, point_b_xyz, point_a_uvw, point_b_uvw, p_u, p_v, num_knot_span_u, num_knot_span_v);
}
else if( local_space_dimension == 3) {
SizeType p_u = mParameters["polynomial_order"].GetArrayItem(0).GetInt();
Expand All @@ -77,7 +94,7 @@ namespace Kratos
SizeType num_knot_span_w = mParameters["number_of_knot_spans"].GetArrayItem(2).GetInt();
SizeType num_knot_span_v = mParameters["number_of_knot_spans"].GetArrayItem(1).GetInt();

CreateAndAddRegularGrid3D(model_part, point_a, point_b, p_u, p_v, p_w, num_knot_span_u, num_knot_span_v, num_knot_span_w);
CreateAndAddRegularGrid3D(model_part, point_a_xyz, point_b_xyz, point_a_uvw, point_b_uvw, p_u, p_v, p_w, num_knot_span_u, num_knot_span_v, num_knot_span_w);
}
else {
KRATOS_ERROR << "Nurbs Geometry Modeler is only available for surfaces and volumes." << std::endl;
Expand All @@ -88,26 +105,30 @@ namespace Kratos
///@}
///@name Private Operations
///@{
void NurbsGeometryModeler::CreateAndAddRegularGrid2D( ModelPart& r_model_part, const Point& A, const Point& B,
SizeType OrderU, SizeType OrderV,SizeType NumKnotSpansU, SizeType NumKnotSpansV)
void NurbsGeometryModeler::CreateAndAddRegularGrid2D( ModelPart& r_model_part, const Point& A_xyz, const Point& B_xyz,
const Point& A_uvw, const Point& B_uvw, SizeType OrderU, SizeType OrderV,SizeType NumKnotSpansU, SizeType NumKnotSpansV)
{
KRATOS_ERROR_IF( B.X() <= A.X() || B.Y() <= A.Y() ) << "NurbsGeometryModeler: "
<< "The two Points A and B must meet the following requirement: (B-A) > (0,0,0). However, (B-A)=" << B-A << std::endl;
KRATOS_ERROR_IF( B_xyz.X() <= A_xyz.X() || B_xyz.Y() <= A_xyz.Y() ) << "NurbsGeometryModeler: "
<< "The two Points A_xyz and B_xyz must meet the following requirement: (B_xyz-A_xyz) > (0,0,0). However, (B_xyz-A_xyz)=" << B_xyz-A_xyz << std::endl;

KRATOS_ERROR_IF( B_uvw.X() <= A_uvw.X() || B_uvw.Y() <= A_uvw.Y() ) << "NurbsGeometryModeler: "
<< "The two Points A_uvw and B_uvw must meet the following requirement: (B_uvw-A_uvw) > (0,0,0). However, (B_uvw-A_uvw)=" << B_uvw-A_uvw << std::endl;


PointerVector<NodeType> points;
double delta_x = std::abs(B.X()-A.X())/OrderU;
double delta_y = std::abs(B.Y()-A.Y())/OrderV;
double delta_x = std::abs(B_xyz.X()-A_xyz.X())/OrderU;
double delta_y = std::abs(B_xyz.Y()-A_xyz.Y())/OrderV;

double x = A.X();
double y = A.Y();
double x = A_xyz.X();
double y = A_xyz.Y();

for( IndexType j=0; j < OrderV+1; ++j){
for( IndexType i = 0; i < OrderU+1; ++i){
// Proper node ID is defined later.
points.push_back(Kratos::make_intrusive<NodeType>(0, x, y, 0.0 ));
x += delta_x;
}
x = A.X();
x = A_xyz.X();
y += delta_y;
}

Expand All @@ -116,18 +137,18 @@ namespace Kratos
Vector knot_vector_u(number_of_knots_u);
for( IndexType i=0; i < number_of_knots_u; ++i){
if( i < OrderU ){
knot_vector_u[i] = 0.0;
knot_vector_u[i] = A_uvw[0];
} else {
knot_vector_u[i] = 1.0;
knot_vector_u[i] = B_uvw[0];
}
}
SizeType number_of_knots_v = 2*OrderV;
Vector knot_vector_v(number_of_knots_v);
for( IndexType i=0; i < number_of_knots_v; ++i){
if( i < OrderV ){
knot_vector_v[i] = 0.0;
knot_vector_v[i] = A_uvw[1];
} else {
knot_vector_v[i] = 1.0;
knot_vector_v[i] = B_uvw[1];
}
}

Expand All @@ -136,16 +157,16 @@ namespace Kratos
points, OrderU, OrderV, knot_vector_u, knot_vector_v);

// Set up knots for knot refinement according to the given number of elements in each direction.
double delta_knot_u = 1.0 / NumKnotSpansU;
double knot_u = 0.0;
double delta_knot_u = (B_uvw[0]-A_uvw[0]) / NumKnotSpansU;
double knot_u = A_uvw[0];
std::vector<double> insert_knots_u(NumKnotSpansU-1);
for( IndexType i = 0; i < NumKnotSpansU-1; i++){
knot_u += delta_knot_u;
insert_knots_u[i] = knot_u;
}

double delta_knot_v = 1.0 / NumKnotSpansV;
double knot_v = 0.0;
double delta_knot_v = (B_uvw[1]-A_uvw[1]) / NumKnotSpansV;
double knot_v = A_uvw[1];
std::vector<double> insert_knots_v(NumKnotSpansV-1);
for( IndexType i = 0; i < NumKnotSpansV-1; i++){
knot_v += delta_knot_v;
Expand Down Expand Up @@ -219,60 +240,64 @@ namespace Kratos
}


void NurbsGeometryModeler::CreateAndAddRegularGrid3D( ModelPart& r_model_part, const Point& A, const Point& B, SizeType OrderU, SizeType OrderV, SizeType OrderW,
void NurbsGeometryModeler::CreateAndAddRegularGrid3D( ModelPart& r_model_part, const Point& A_xyz, const Point& B_xyz,
const Point& A_uvw, const Point& B_uvw, SizeType OrderU, SizeType OrderV, SizeType OrderW,
SizeType NumKnotSpansU, SizeType NumKnotSpansV, SizeType NumKnotSpansW )
{
KRATOS_ERROR_IF( B.X() <= A.X() || B.Y() <= A.Y() || B.Z() <= A.Z() ) << "NurbsGeometryModeler: "
<< "The two Points A and B must meet the following requirement: (B-A) > (0,0,0). However, (B-A)=" << B-A << std::endl;
KRATOS_ERROR_IF( B_xyz.X() <= A_xyz.X() || B_xyz.Y() <= A_xyz.Y() || B_xyz.Z() <= A_xyz.Z() ) << "NurbsGeometryModeler: "
<< "The two Points A_xyz and B_xyz must meet the following requirement: (B_xyz-A_xyz) > (0,0,0). However, (B_xyz-A_xyz)=" << B_xyz-A_xyz << std::endl;

KRATOS_ERROR_IF( B_uvw.X() <= A_uvw.X() || B_uvw.Y() <= A_uvw.Y() || B_uvw.Z() <= A_uvw.Z() ) << "NurbsGeometryModeler: "
<< "The two Points A_uvw and B_uvw must meet the following requirement: (B_uvw-A_uvw) > (0,0,0). However, (B_uvw-A_uvw)=" << B_uvw-A_uvw << std::endl;

// Set up control points.
// Note: The CP's are uniformly subdivided with increasing p.
// This is allowed, as all lines of the cartesian grid are straight.
PointerVector<NodeType> points;
double delta_x = std::abs(B.X()-A.X())/OrderU;
double delta_y = std::abs(B.Y()-A.Y())/OrderV;
double delta_z = std::abs(B.Z()-A.Z())/OrderW;
double x = A.X();
double y = A.Y();
double z = A.Z();
double delta_x = std::abs(B_xyz.X()-A_xyz.X())/OrderU;
double delta_y = std::abs(B_xyz.Y()-A_xyz.Y())/OrderV;
double delta_z = std::abs(B_xyz.Z()-A_xyz.Z())/OrderW;
double x = A_xyz.X();
double y = A_xyz.Y();
double z = A_xyz.Z();
for( IndexType k=0; k < OrderW+1; ++k){
for( IndexType j=0; j < OrderV+1; ++j){
for( IndexType i = 0; i < OrderU+1; ++i){
points.push_back(Kratos::make_intrusive<NodeType>(0, x, y, z ));
x += delta_x;
}
x = A.X();
x = A_xyz.X();
y += delta_y;
}
y = A.Y();
y = A_xyz.Y();
z += delta_z;
}
// Set up knot vectors according to the given polynomial degree p.
SizeType number_of_knots_u = 2*OrderU;
Vector knot_vector_u(number_of_knots_u);
for( IndexType i=0; i < number_of_knots_u; ++i){
if( i < OrderU ){
knot_vector_u[i] = 0.0;
knot_vector_u[i] = A_uvw[0];
} else {
knot_vector_u[i] = 1.0;
knot_vector_u[i] = B_uvw[0];
}
}
SizeType number_of_knots_v = 2*OrderV;
Vector knot_vector_v(number_of_knots_v);
for( IndexType i=0; i < number_of_knots_v; ++i){
if( i < OrderV ){
knot_vector_v[i] = 0.0;
knot_vector_v[i] = A_uvw[1];
} else {
knot_vector_v[i] = 1.0;
knot_vector_v[i] = B_uvw[1];
}
}
SizeType number_of_knots_w = 2*OrderW;
Vector knot_vector_w(number_of_knots_w);
for( IndexType i=0; i < number_of_knots_w; ++i){
if( i < OrderW ){
knot_vector_w[i] = 0.0;
knot_vector_w[i] = A_uvw[2];
} else {
knot_vector_w[i] = 1.0;
knot_vector_w[i] = B_uvw[2];
}
}

Expand All @@ -281,24 +306,24 @@ namespace Kratos
points, OrderU, OrderV, OrderW, knot_vector_u, knot_vector_v, knot_vector_w);

// Set up knots for knot refinement according to the given number of elements in each direction.
double delta_knot_u = 1.0 / NumKnotSpansU;
double knot_u = 0.0;
double delta_knot_u = (B_uvw[0]-A_uvw[0]) / NumKnotSpansU;
double knot_u = A_uvw[0];
std::vector<double> insert_knots_u(NumKnotSpansU-1);
for( IndexType i = 0; i < NumKnotSpansU-1; i++){
knot_u += delta_knot_u;
insert_knots_u[i] = knot_u;
}

double delta_knot_v = 1.0 / NumKnotSpansV;
double knot_v = 0.0;
double delta_knot_v = (B_uvw[1]-A_uvw[1]) / NumKnotSpansV;
double knot_v = A_uvw[1];
std::vector<double> insert_knots_v(NumKnotSpansV-1);
for( IndexType i = 0; i < NumKnotSpansV-1; i++){
knot_v += delta_knot_v;
insert_knots_v[i] = knot_v;
}

double delta_knot_w = 1.0 / NumKnotSpansW;
double knot_w = 0.0;
double delta_knot_w = (B_uvw[2]-A_uvw[2])/ NumKnotSpansW;
double knot_w = A_uvw[2];
std::vector<double> insert_knots_w(NumKnotSpansW-1);
for( IndexType i = 0; i < NumKnotSpansW-1; i++){
knot_w += delta_knot_w;
Expand Down
8 changes: 4 additions & 4 deletions applications/IgaApplication/custom_modelers/nurbs_geometry_modeler.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -100,8 +100,8 @@ class KRATOS_API(IGA_APPLICATION) NurbsGeometryModeler
* @param NumKnotSpans Number of equidistant elements/knot spans in each direction u,v.
* @note The CP'S are defined as nodes and added to the rModelPart.
**/
void CreateAndAddRegularGrid2D( ModelPart& r_model_part, const Point& A, const Point& B, SizeType OrderU, SizeType OrderV,
SizeType NumKnotSpansU, SizeType NumKnotSpansV );
void CreateAndAddRegularGrid2D( ModelPart& r_model_part, const Point& A_xyz, const Point& B_xyz, const Point& A_uvw, const Point& B_uvw,
SizeType OrderU, SizeType OrderV, SizeType NumKnotSpansU, SizeType NumKnotSpansV );

/**
* @brief Creates a cartesian grid composed out of trivariant B-spline cubes.
Expand All @@ -111,8 +111,8 @@ class KRATOS_API(IGA_APPLICATION) NurbsGeometryModeler
* @param NumKnotSpans Number of equidistant elements/knot spans in each direction u,v,w.
* @note The CP'S are defined as nodes and added to the rModelPart.
**/
void CreateAndAddRegularGrid3D( ModelPart& r_model_part, const Point& A, const Point& B, SizeType OrderU, SizeType OrderV, SizeType OrderW,
SizeType NumKnotSpansU, SizeType NumKnotSpansV, SizeType NumKnotSpansW );
void CreateAndAddRegularGrid3D( ModelPart& r_model_part, const Point& A_xyz, const Point& B_xyz, const Point& A_uvw, const Point& B_uvw,
SizeType OrderU, SizeType OrderV, SizeType OrderW, SizeType NumKnotSpansU, SizeType NumKnotSpansV, SizeType NumKnotSpansW );

};

Expand Down
18 changes: 15 additions & 3 deletions ...Application/custom_processes/assign_integration_points_to_background_elements_process.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@

// Project includes
#include "assign_integration_points_to_background_elements_process.h"
#include "geometries/nurbs_volume_geometry.h"

namespace Kratos
{
Expand Down Expand Up @@ -52,6 +53,14 @@ namespace Kratos
const CoordinatesArrayType lower_point = p_geometry->begin()->GetInitialPosition();
const CoordinatesArrayType upper_point = (p_geometry->end()-1)->GetInitialPosition();

const auto p_nurbs_volume_geo = dynamic_pointer_cast<NurbsVolumeGeometry<PointerVector<Node>>>(p_geometry);
const auto knots_u = p_nurbs_volume_geo->KnotsU();
const SizeType n_knots_u = p_nurbs_volume_geo->NumberOfKnotsU();
const auto knots_v = p_nurbs_volume_geo->KnotsV();
const SizeType n_knots_v = p_nurbs_volume_geo->NumberOfKnotsV();
const auto knots_w = p_nurbs_volume_geo->KnotsW();
const SizeType n_knots_w = p_nurbs_volume_geo->NumberOfKnotsW();

// Create quadrature points in parameter space of nurbs volume
IntegrationPointsArrayType integration_points(embedded_model_part.NumberOfElements());
const auto element_itr_begin = embedded_model_part.ElementsBegin();
Expand All @@ -65,9 +74,12 @@ namespace Kratos

// Map point into parameter space
CoordinatesArrayType local_point;
local_point[0] = (ip_physical_space[0] - lower_point[0]) / std::abs( lower_point[0] - upper_point[0]);
local_point[1] = (ip_physical_space[1] - lower_point[1]) / std::abs( lower_point[1] - upper_point[1]);
local_point[2] = (ip_physical_space[2] - lower_point[2]) / std::abs( lower_point[2] - upper_point[2]);
local_point[0] = ((ip_physical_space[0] - lower_point[0]) / std::abs( lower_point[0] - upper_point[0])
* std::abs(knots_u[n_knots_u-1] - knots_u[0])) + knots_u[0];
local_point[1] = ((ip_physical_space[1] - lower_point[1]) / std::abs( lower_point[1] - upper_point[1])
* std::abs(knots_v[n_knots_v-1] - knots_v[0])) + knots_v[0];
local_point[2] = ((ip_physical_space[2] - lower_point[2]) / std::abs( lower_point[2] - upper_point[2])
* std::abs(knots_w[n_knots_w-1] - knots_w[0])) + knots_w[0];

integration_points[i] = IntegrationPoint<3>(local_point, tmp_integration_points[0].Weight());
}
Expand Down
17 changes: 14 additions & 3 deletions ...IgaApplication/custom_processes/map_nurbs_volume_results_to_embedded_geometry_process.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,14 +80,25 @@ namespace Kratos
const CoordinatesArrayType lower_point = p_geometry->begin()->GetInitialPosition();
const CoordinatesArrayType upper_point = (p_geometry->end()-1)->GetInitialPosition();

const auto p_nurbs_volume_geo = dynamic_pointer_cast<NurbsVolumeGeometry<PointerVector<Node>>>(p_geometry);
const auto knots_u = p_nurbs_volume_geo->KnotsU();
const SizeType n_knots_u = p_nurbs_volume_geo->NumberOfKnotsU();
const auto knots_v = p_nurbs_volume_geo->KnotsV();
const SizeType n_knots_v = p_nurbs_volume_geo->NumberOfKnotsV();
const auto knots_w = p_nurbs_volume_geo->KnotsW();
const SizeType n_knots_w = p_nurbs_volume_geo->NumberOfKnotsW();

const auto node_itr_begin = embedded_model_part.NodesBegin();
IndexPartition<IndexType>(embedded_model_part.NumberOfNodes()).for_each([&](IndexType i) {
auto node_itr = node_itr_begin + i;
// Map point into parameter space
CoordinatesArrayType local_point;
local_point[0] = (node_itr->X() - lower_point[0]) / std::abs( lower_point[0] - upper_point[0]);
local_point[1] = (node_itr->Y() - lower_point[1]) / std::abs( lower_point[1] - upper_point[1]);
local_point[2] = (node_itr->Z() - lower_point[2]) / std::abs( lower_point[2] - upper_point[2]);
local_point[0] = ((node_itr->X() - lower_point[0]) / std::abs( lower_point[0] - upper_point[0])
* std::abs(knots_u[n_knots_u-1] - knots_u[0])) + knots_u[0];
local_point[1] = ((node_itr->Y() - lower_point[1]) / std::abs( lower_point[1] - upper_point[1])
* std::abs(knots_v[n_knots_v-1] - knots_v[0])) + knots_v[0];
local_point[2] = ((node_itr->Z() - lower_point[2]) / std::abs( lower_point[2] - upper_point[2])
* std::abs(knots_w[n_knots_w-1] - knots_w[0])) + knots_w[0];
integration_points[i] = IntegrationPoint<3>(local_point, 1.0);
});

Expand Down
6 changes: 4 additions & 2 deletions ...ations/IgaApplication/tests/test_map_nurbs_volume_results_to_embedded_geometry_process.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,8 +30,10 @@ def test_MapNurbsVolumeResultsToEmbeddedGeometryProcess(self):
"Parameters": {
"model_part_name" : "NurbsMesh",
"geometry_name" : "NurbsVolume",
"lower_point": [-2.1, -2.1, -2.1],
"upper_point": [2.1, 2.1, 12.1],
"lower_point_xyz": [-2.1, -2.1, -2.1],
"upper_point_xyz": [2.1, 2.1, 12.1],
"lower_point_uvw": [0.0, 0.0, 0.0],
"upper_point_uvw": [1.0, 1.0, 1.0],
"polynomial_order" : [1, 2, 2],
"number_of_knot_spans" : [1,2,2]}
},{
Expand Down
Loading

0 comments on commit ca793fd

Please sign in to comment.