Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Feature curved pyramids #1036

Open
wants to merge 21 commits into
base: main
Choose a base branch
from
Open

Feature curved pyramids #1036

Show file tree
Hide file tree
Changes from 20 commits
Commits
Show all changes
21 commits
Select commit Hold shift + click to select a range
bd3420e
added numbering for pyramid element in t8_eclass.c
jfussbro Apr 5, 2024
f40381a
added cad pyramid test
jfussbro Apr 5, 2024
34f4664
added comment with link to issue
jfussbro Apr 5, 2024
ea35996
Merge remote-tracking branch 'origin/main' into feature_curved_pyramids
jfussbro Apr 9, 2024
2bdcf9d
Merge branch 'feature-puma_patches_batch_processing_geometry3' into t…
jfussbro Apr 10, 2024
10f5526
Merge remote-tracking branch 'origin/main' into feature_curved_pyramids
jfussbro Apr 15, 2024
16ff7e1
Corrected fault in lookup table
jfussbro Apr 15, 2024
01cdd29
Added mapping algorithm for cad pyramids
jfussbro Apr 15, 2024
aaa3911
Adapted test for cad pyramids to now existing mapping algorithm
jfussbro Apr 15, 2024
108b2c1
Merge branch 'feature-puma_patches_batch_processing_geometry3' into f…
jfussbro Apr 15, 2024
3e493c6
Merge branch 'feature-puma_patches_batch_processing_geometry3' into f…
jfussbro Apr 23, 2024
5f71d4d
enabled read-in of cad-linked pyramids
jfussbro Apr 23, 2024
0dae72b
bug fix - wrong scaling factor of edges of pyramid face 4
jfussbro Apr 24, 2024
156f53c
added element centroid as testing point
jfussbro Apr 24, 2024
62bf412
Merge remote-tracking branch 'origin/main' into feature_curved_pyramids
jfussbro Apr 29, 2024
a95045f
bug fix
jfussbro Apr 29, 2024
9af614c
Merge remote-tracking branch 'origin/main' into feature_curved_pyramids
jfussbro May 22, 2024
4b69d74
indent
jfussbro May 22, 2024
a10bae6
removed old todo
jfussbro May 22, 2024
cbfd3f4
indent
jfussbro May 22, 2024
1d645c9
added suggestions from code review
jfussbro May 24, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 5 additions & 5 deletions src/t8_cmesh/t8_cmesh_readmshfile.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1016,11 +1016,11 @@ t8_cmesh_msh_file_4_read_eles (t8_cmesh_t cmesh, FILE *fp, sc_hash_t *vertices,
T8_ASSERT (cad_geometry_base->t8_geom_get_type () == T8_GEOMETRY_TYPE_CAD);
const t8_geometry_cad_c *cad_geometry = dynamic_cast<const t8_geometry_cad_c *> (cad_geometry_base);
/* Check for right element class */
if (eclass != T8_ECLASS_TRIANGLE && eclass != T8_ECLASS_QUAD && eclass != T8_ECLASS_HEX
&& eclass != T8_ECLASS_TET && eclass != T8_ECLASS_PRISM) {
t8_errorf (
"%s element detected. The cad geometry currently only supports quad, tri, hex and prism elements.",
t8_eclass_to_string[eclass]);
if (eclass != T8_ECLASS_TRIANGLE && eclass != T8_ECLASS_QUAD && eclass != T8_ECLASS_TET
&& eclass != T8_ECLASS_HEX && eclass != T8_ECLASS_PRISM && eclass != T8_ECLASS_PYRAMID) {
t8_errorf ("%s element detected. The occ geometry currently only supports quad, tri, tet, hex, prism and "
"pyramid elements.",
t8_eclass_to_string[eclass]);
goto die_ele;
}
int tree_is_linked = 0;
Expand Down
15 changes: 6 additions & 9 deletions src/t8_eclass.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,6 @@ const int t8_face_vertex_to_tree_vertex[T8_ECLASS_COUNT][T8_ECLASS_MAX_FACES][T8
{ { 0, 2, 4 }, { 1, 3, 4 }, { 0, 1, 4 }, { 2, 3, 4 }, { 0, 1, 2, 3 } } /* pyramid */
};

/* TODO: tet, pyramid */
const int t8_face_edge_to_tree_edge[T8_ECLASS_COUNT][T8_ECLASS_MAX_FACES][T8_ECLASS_MAX_EDGES_2D] = {
{ { -1 } }, /* vertex */
{ { 0 } }, /* line */
Expand All @@ -48,22 +47,21 @@ const int t8_face_edge_to_tree_edge[T8_ECLASS_COUNT][T8_ECLASS_MAX_FACES][T8_ECL
{ { 8, 10, 4, 6 }, { 9, 11, 5, 7 }, { 8, 9, 0, 2 }, { 10, 11, 1, 3 }, { 4, 5, 0, 1 }, { 6, 7, 2, 3 } }, /* hex */
{ { 3, 4, 5 }, { 1, 2, 5 }, { 0, 2, 4 }, { 0, 1, 3 } }, /* tet */
{ { 0, 7, 3, 6 }, { 1, 8, 4, 7 }, { 2, 6, 5, 8 }, { 0, 1, 2 }, { 3, 4, 5 } }, /* prism */
{ { -1 } }, /* pyramid */
{ { 0, 4, 7 }, { 1, 5, 6 }, { 2, 4, 5 }, { 3, 6, 7 }, { 0, 1, 2, 3 } } /* pyramid */
};

/* TODO: tet, pyramid */
const int t8_face_to_edge_neighbor[T8_ECLASS_COUNT][T8_ECLASS_MAX_FACES][T8_ECLASS_MAX_CORNERS_2D] = {
{ { -1 } }, /* vertex */
{ { -1 } }, /* line */
{ { 2, 3 }, { 2, 3 }, { 0, 1 }, { 0, 1 } }, /* quad */
{ { 2, 1 }, { 2, 0 }, { 1, 0 } }, /* triangle */
{ { 0, 1, 2, 3 }, { 0, 1, 2, 3 }, { 4, 5, 6, 7 }, { 4, 5, 6, 7 }, { 8, 9, 10, 11 }, { 8, 9, 10, 11 } }, /* hex */
{ { 0, 1, 2 }, { 0, 3, 4 }, { 1, 3, 5 }, { 2, 4, 5 } }, /* tet */
{ { 1, 2, 4, 5 }, { 0, 2, 3, 5 }, { 0, 1, 3, 4 }, { 6, 7, 8 }, { 6, 7, 8 } }, /* prism */
{ { -1 } }, /* pyramid */
{ { 1, 2, 4, 5 }, { 0, 2, 3, 5 }, { 0, 1, 3, 4 }, { 6, 7, 8 }, { 6, 7, 8 } },
/* prism */ /* prism */
{ { 2, 3, 5, 6 }, { 2, 3, 4, 7 }, { 0, 1, 6, 7 }, { 0, 1, 4, 5 }, { 4, 5, 6, 7 } } /* pyramid */
};

/* TODO: tet, pyramid */
const int t8_edge_vertex_to_tree_vertex[T8_ECLASS_COUNT][T8_ECLASS_MAX_EDGES][2] = {
{ { -1 } }, /* vertex */
{ { 0 }, { 1 } }, /* line */
Expand All @@ -83,10 +81,9 @@ const int t8_edge_vertex_to_tree_vertex[T8_ECLASS_COUNT][T8_ECLASS_MAX_EDGES][2]
{ 3, 7 } }, /* hex */
{ { 0, 1 }, { 0, 2 }, { 0, 3 }, { 1, 2 }, { 1, 3 }, { 2, 3 } }, /* tet */
{ { 1, 2 }, { 0, 2 }, { 0, 1 }, { 4, 5 }, { 3, 5 }, { 3, 4 }, { 1, 4 }, { 2, 5 }, { 0, 3 } }, /* prism */
{ { -1 } }, /* pyramid */
{ { 0, 2 }, { 1, 3 }, { 0, 1 }, { 2, 3 }, { 0, 4 }, { 1, 4 }, { 3, 4 }, { 2, 4 } } /* pyramid */
};

/* TODO: tet, pyramid */
const int t8_edge_to_face[T8_ECLASS_COUNT][T8_ECLASS_MAX_EDGES][2] = {
{ { -1 } }, /* vertex */
{ { 0 } }, /* line */
Expand All @@ -106,7 +103,7 @@ const int t8_edge_to_face[T8_ECLASS_COUNT][T8_ECLASS_MAX_EDGES][2] = {
{ 1, 3 } }, /* hex */
{ { 2, 3 }, { 1, 3 }, { 1, 2 }, { 0, 3 }, { 0, 2 }, { 0, 1 } }, /* tet */
{ { 0, 3 }, { 1, 3 }, { 2, 3 }, { 0, 4 }, { 1, 4 }, { 2, 4 }, { 0, 2 }, { 0, 1 }, { 1, 2 } }, /* prism */
{ { -1 } }, /* pyramid */
{ { 0, 4 }, { 1, 4 }, { 2, 4 }, { 3, 4 }, { 0, 2 }, { 1, 2 }, { 1, 3 }, { 0, 3 } } /* pyramid */
};

const int t8_eclass_face_orientation[T8_ECLASS_COUNT][T8_ECLASS_MAX_FACES] = {
Expand Down
71 changes: 71 additions & 0 deletions src/t8_geometry/t8_geometry_helpers.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -411,6 +411,77 @@ t8_geom_get_scaling_factor_of_edge_on_face_tet (const int edge, const int face,
return (1.0 - (orthogonal_vector[edge][face] / max_orthogonal_vector[edge][face]));
}

double
t8_geom_get_scaling_factor_of_edge_on_face_pyramid (const int edge, const int face, const double *ref_coords)
{
/* Save the orthogonal direction and the maximum of that direction
* of an edge in reference space on one of the neighbouring faces.
* /|
* / |
* / |
* / |
* / |--edge
* / --|--face
* / |
* / <~~| orthogonal direction
* /<----o--| maximum othogonal direction
* /_________|
*/
const double orthogonal_direction[8][5]
= { { ref_coords[2], 0, 0, 0, ref_coords[0] },
{ 0, ref_coords[2], 0, 0, (1 - ref_coords[0]) },
{ 0, 0, ref_coords[2], 0, ref_coords[1] },
{ 0, 0, 0, ref_coords[2], (1 - ref_coords[1]) },
{ (ref_coords[1] - ref_coords[2]), 0, (ref_coords[0] - ref_coords[2]), 0, 0 },
{ 0, (ref_coords[1] - ref_coords[2]), (1 - ref_coords[0]), 0, 0 },
{ 0, (1 - ref_coords[1]), 0, (1 - ref_coords[0]), 0 },
{ (1 - ref_coords[1]), 0, 0, (ref_coords[0] - ref_coords[2]), 0 } };
const double max_orthogonal_direction[8][5] = { { ref_coords[1], 0, 0, 0, 1 },
{ 0, ref_coords[1], 0, 0, 1 },
{ 0, 0, ref_coords[0], 0, 1 },
{ 0, 0, 0, ref_coords[0], 1 },
{ (1 - ref_coords[2]), 0, (1 - ref_coords[0]), 0, 0 },
{ 0, (1 - ref_coords[2]), (1 - ref_coords[2]), 0, 0 },
{ 0, (1 - ref_coords[2]), 0, (1 - ref_coords[2]), 0 },
{ (1 - ref_coords[2]), 0, 0, (1 - ref_coords[2]), 0 } };
Comment on lines +430 to +446
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Most of the values are calculated and not used

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Every value which is not 0 will get used at some point. Removing the 0's is not possible, since calling the lookup with the edge number wouldn't be possible anymore


/* Check that the combination of edge and face is valid */
T8_ASSERT (face == t8_edge_to_face[T8_ECLASS_PYRAMID][edge][0]
|| face == t8_edge_to_face[T8_ECLASS_PYRAMID][edge][1]);

/* If the maximum orthogonal direction is 1, the reference coordinate lies on
* one of the edge nodes and the scaling factor is therefore 0, because the displacement
* at the nodes is always 0.
* In all other cases the scaling factor is determined with one minus the relation of the orthogonal direction
* to the maximum orthogonal direction. */
if (max_orthogonal_direction[edge][face] == 0) {
return 0;
}
else {
return (1.0 - (orthogonal_direction[edge][face] / max_orthogonal_direction[edge][face]));
}
}

double
t8_geom_get_scaling_factor_face_through_volume_pyramid (const int face, const double *ref_coords)
{
/* The function computes the scaling factor of any displacement of a pyramid face, throughout the
* volume of the element. The scaling factor is calculated accordingly to
* t8_geom_get_scaling_factor_of_edge_on_face_pyramid with 1 - (orthogonal_direction / max_orthogonal_direction) */

const double orthogonal_direction[5] = { (ref_coords[0] - ref_coords[2]), (1 - ref_coords[0]),
(ref_coords[1] - ref_coords[2]), (1 - ref_coords[1]), ref_coords[2] };
const double max_orthogonal_direction[5]
= { (1 - ref_coords[2]), (1 - ref_coords[2]), (1 - ref_coords[2]), (1 - ref_coords[2]),
(ref_coords[0] >= ref_coords[1] ? ref_coords[1] : ref_coords[0]) };
Comment on lines +472 to +476
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same here

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why are they not used?


if (max_orthogonal_direction[face] == 0) {
return 1.0;
}

return (1.0 - (orthogonal_direction[face] / max_orthogonal_direction[face]));
}

void
t8_geom_get_tet_face_intersection (const int face, const double *ref_coords, double face_intersection[3])
{
Expand Down
19 changes: 19 additions & 0 deletions src/t8_geometry/t8_geometry_helpers.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -163,6 +163,25 @@ t8_geom_get_tet_face_intersection (const int face_index, const double *ref_coord
double
t8_geom_get_scaling_factor_of_edge_on_face_prism (int edge_index, int face_index, const double *ref_coords);

/** Calculates the scaling factor for the displacement of an edge over a face of a pyramid element.
* \param [in] edge_index Index of the edge, whose displacement should be scaled.
* \param [in] face_index Index of the face, the displacement should be scaled on.
* \param [in] ref_coords Array containing the coordinates of the reference point.
* \return The scaling factor of the edge displacement on the face
* at the point of the reference coordinates.
*/
double
t8_geom_get_scaling_factor_of_edge_on_face_pyramid (int edge_index, int face_index, const double *ref_coords);

/** Calculates the scaling factor for the displacement of an face through the volume of a pyramid element.
* \param [in] face_index Index of the displaced face.
* \param [in] ref_coords Array containing the coordinates of the reference point.
* \return The scaling factor of the face displacement
* at the point of the reference coordinates inside the pyramid volume.
*/
double
t8_geom_get_scaling_factor_face_through_volume_pyramid (const int face, const double *ref_coords);

/** Calculates the scaling factor for the displacement of an face through the volume of a prism element.
* \param [in] face_index Index of the displaced face.
* \param [in] ref_coords Array containing the coordinates of the reference point.
Expand Down
Loading