Skip to content

GitLab

Commit 6e849737 authored by Daniel Drzisga's avatar Daniel Drzisga
Browse files

finished p2 pull_edges, p2 pull_vertices and p2 interpolation

parent a0eed49b
Hide whitespace changes
Inline Side-by-side
include/p2functionspace/p2edge.hpp
View file @ 6e849737
......@@ -34,14 +34,14 @@ inline void free(Edge& edge, size_t memory_id, size_t minLevel, size_t maxLevel)
}
}
inline void interpolate(Edge& edge, size_t memory_id, std::function<double(const hhg::Point3D&)>& expr, size_t level, bool boundaryDofs)
inline void interpolate(Edge& edge, size_t memory_id, std::function<double(const hhg::Point3D&)>& expr, size_t level, size_t flag)
{
size_t rowsize = levelinfo::num_microvertices_per_edge(level) + levelinfo::num_microedges_per_edge(level);
Point3D x = edge.v0->coords;
Point3D dx = edge.direction / (double) (rowsize - 1);
x += dx;
if (boundaryDofs)
if (testFlag(edge.type, flag))
{
for (size_t i = 1; i < rowsize-1; ++i)
{
......@@ -50,35 +50,37 @@ inline void interpolate(Edge& edge, size_t memory_id, std::function<double(const
}
}
Point3D v(edge.faces[0]->get_vertex_opposite_to_edge(edge)->coords);
Point3D dirv((v - edge.v0->coords) / ((double) rowsize - 1));
x = edge.v0->coords + dx + dirv;
size_t offset = rowsize;
for (size_t i = 1; i < rowsize-1-1; ++i)
{
edge.data[memory_id][level-2][offset+i] = expr(x);
x += dx;
}
offset += rowsize-1 + rowsize-2;
if (edge.faces.size() == 2)
if (testFlag(hhg::Inner, flag))
{
v = edge.faces[1]->get_vertex_opposite_to_edge(edge)->coords;
dirv = (v - edge.v0->coords) / ((double) rowsize - 1);
Point3D v(edge.faces[0]->get_vertex_opposite_to_edge(edge)->coords);
Point3D dirv((v - edge.v0->coords) / ((double) rowsize - 1));
x = edge.v0->coords + dx + dirv;
size_t offset = rowsize;
for (size_t i = 1; i < rowsize-1-1; ++i)
{
edge.data[memory_id][level-2][offset+i] = expr(x);
x += dx;
}
}
offset += rowsize-1 + rowsize-2;
if (edge.faces.size() == 2)
{
v = edge.faces[1]->get_vertex_opposite_to_edge(edge)->coords;
dirv = (v - edge.v0->coords) / ((double) rowsize - 1);
x = edge.v0->coords + dx + dirv;
for (size_t i = 1; i < rowsize-1-1; ++i)
{
edge.data[memory_id][level-2][offset+i] = expr(x);
x += dx;
}
}
}
}
inline void pull_vertices(Edge& edge, size_t memory_id, size_t level)
......
include/p2functionspace/p2face.hpp
View file @ 6e849737
......@@ -53,7 +53,7 @@ inline void interpolate(Face& face, size_t memory_id, std::function<double(const
for (size_t i = 0; i < rowsize-5; ++i)
{
x = x0;
x += (i+1) * d2 *2 + d0 * 2;
x += (i+2) * d2 + d0 * 2;
for (size_t j = 0; j < inner_rowsize-5; ++j)
{
......@@ -84,13 +84,22 @@ inline void print(Face & face, size_t memory_id, size_t level) {
inline void pull_edges(Face& face, size_t memory_id, size_t level)
{
size_t rowsize = levelinfo::num_microvertices_per_edge(level) + levelinfo::num_microedges_per_edge(level);
size_t off_edge_size = rowsize - 1 + rowsize - 2;
double* edge_data_0 = NULL;
double* edge_data_1 = NULL;
// double* edge_data_2 = NULL;
double* edge_data_2 = NULL;
size_t edge_0_size = rowsize + face.edges[0]->faces.size() * off_edge_size;
size_t edge_1_size = rowsize + face.edges[1]->faces.size() * off_edge_size;
size_t edge_2_size = rowsize + face.edges[2]->faces.size() * off_edge_size;
size_t edge_0_offset = rowsize + face.edges[0]->face_index(face) * off_edge_size;
size_t edge_1_offset = rowsize + face.edges[1]->face_index(face) * off_edge_size;
size_t edge_2_offset = rowsize + face.edges[2]->face_index(face) * off_edge_size;
MPI_Request req0;
MPI_Request req1;
// MPI_Request req2;
MPI_Request req2;
int rk = hhg::Comm::get().rk;
......@@ -102,13 +111,13 @@ inline void pull_edges(Face& face, size_t memory_id, size_t level)
}
else
{
MPI_Send(&face.edges[0]->data[memory_id][level-2][0], rowsize, MPI_DOUBLE, face.rank, face.edges[0]->id, MPI_COMM_WORLD);
MPI_Send(&face.edges[0]->data[memory_id][level-2][0], edge_0_size, MPI_DOUBLE, face.rank, face.edges[0]->id, MPI_COMM_WORLD);
}
}
else if (face.rank == rk)
{
edge_data_0 = new double[rowsize];
MPI_Irecv(edge_data_0, rowsize, MPI_DOUBLE, face.edges[0]->rank, face.edges[0]->id, MPI_COMM_WORLD, &req0);
edge_data_0 = new double[edge_0_size];
MPI_Irecv(edge_data_0, edge_0_size, MPI_DOUBLE, face.edges[0]->rank, face.edges[0]->id, MPI_COMM_WORLD, &req0);
}
if (face.edges[1]->rank == rk)
......@@ -119,31 +128,31 @@ inline void pull_edges(Face& face, size_t memory_id, size_t level)
}
else
{
MPI_Send(&face.edges[1]->data[memory_id][level-2][0], rowsize, MPI_DOUBLE, face.rank, face.edges[1]->id, MPI_COMM_WORLD);
MPI_Send(&face.edges[1]->data[memory_id][level-2][0], edge_1_size, MPI_DOUBLE, face.rank, face.edges[1]->id, MPI_COMM_WORLD);
}
}
else if (face.rank == rk)
{
edge_data_1 = new double[rowsize];
MPI_Irecv(edge_data_1, rowsize, MPI_DOUBLE, face.edges[1]->rank, face.edges[1]->id, MPI_COMM_WORLD, &req1);
edge_data_1 = new double[edge_1_size];
MPI_Irecv(edge_data_1, edge_1_size, MPI_DOUBLE, face.edges[1]->rank, face.edges[1]->id, MPI_COMM_WORLD, &req1);
}
// if (face.edges[2]->rank == rk)
// {
// if (face.rank == rk)
// {
// edge_data_2 = face.edges[2]->data[memory_id][level-2];
// }
// else
// {
// MPI_Send(&face.edges[2]->data[memory_id][level-2][0], rowsize, MPI_DOUBLE, face.rank, face.edges[2]->id, MPI_COMM_WORLD);
// }
// }
// else if (face.rank == rk)
// {
// edge_data_2 = new double[rowsize];
// MPI_Irecv(edge_data_2, rowsize, MPI_DOUBLE, face.edges[2]->rank, face.edges[2]->id, MPI_COMM_WORLD, &req2);
// }
if (face.edges[2]->rank == rk)
{
if (face.rank == rk)
{
edge_data_2 = face.edges[2]->data[memory_id][level-2];
}
else
{
MPI_Send(&face.edges[2]->data[memory_id][level-2][0], edge_2_size, MPI_DOUBLE, face.rank, face.edges[2]->id, MPI_COMM_WORLD);
}
}
else if (face.rank == rk)
{
edge_data_2 = new double[edge_2_size];
MPI_Irecv(edge_data_2, edge_2_size, MPI_DOUBLE, face.edges[2]->rank, face.edges[2]->id, MPI_COMM_WORLD, &req2);
}
if (face.rank == rk)
{
......@@ -159,29 +168,38 @@ inline void pull_edges(Face& face, size_t memory_id, size_t level)
MPI_Wait(&req1, MPI_STATUS_IGNORE);
}
// if (face.edges[2]->rank != rk)
// {
// MPI_Wait(&req2, MPI_STATUS_IGNORE);
// }
if (face.edges[2]->rank != rk)
{
MPI_Wait(&req2, MPI_STATUS_IGNORE);
}
// edge 0
if (face.edge_orientation[0] == 1)
{
for (size_t i = 0; i < rowsize + rowsize - 1; ++i)
for (size_t i = 0; i < rowsize; ++i)
{
// fmt::print("idx[{}] = {}\n", i, i);
face_data[i] = edge_data_0[i];
}
for (size_t i = 1; i < rowsize-2; ++i)
{
// fmt::print("idx[{}] = {}\n", i, rowsize + i);
face_data[rowsize + i] = edge_data_0[edge_0_offset + i];
}
}
else
{
for (size_t i = 0; i < rowsize; ++i)
{
// fmt::print("idx[{}] = {}\n", i, rowsize - 1 - i);
face_data[i] = edge_data_0[rowsize - 1 - i];
}
for (size_t i = 0; i < rowsize-1; ++i)
for (size_t i = 1; i < rowsize-2; ++i)
{
face_data[i + rowsize] = edge_data_0[(rowsize + rowsize - 1) - 1 - i];
// fmt::print("idx[{}] = {}\n", i, (rowsize - 1) - 1 - i);
face_data[rowsize + i] = edge_data_0[edge_0_offset + (rowsize - 1) - 1 - i];
}
}
......@@ -196,59 +214,86 @@ inline void pull_edges(Face& face, size_t memory_id, size_t level)
size_t idx = rowsize - 1;
for (size_t i = 0; i < rowsize; ++i)
{
// fmt::print("idx[{}] = {}\n", i, idx);
face_data[idx] = edge_data_1[i];
idx += rowsize - 1 - i;
fmt::print("idx = {}\n", idx);
}
idx = rowsize + rowsize - 1 - 2;
for (size_t i = 1; i < rowsize-2; ++i)
{
// fmt::print("idx[{}] = {}\n", i, idx);
face_data[idx] = edge_data_1[edge_1_offset + i];
idx += rowsize - 1 - i;
}
}
else
{
size_t idx = levelinfo::num_microvertices_per_face(level) + levelinfo::num_microedges_per_face(level) - 1;
for (size_t i = 0; i < rowsize; ++i)
{
// fmt::print("idx[{}] = {}\n", i, idx);
face_data[idx] = edge_data_1[i];
idx += rowsize - 1 - i;
idx -= i + 1;
}
fmt::print("idx2 = {}\n", idx);
idx = levelinfo::num_microvertices_per_face(level) + levelinfo::num_microedges_per_face(level) - 1 - 4;
for (size_t i = 1; i < rowsize-2; ++i)
{
// fmt::print("idx[{}] = {}\n", i, idx);
face_data[idx] = edge_data_1[edge_1_offset + i];
idx -= i + 2;
}
}
if (face.edges[1]->rank != rk)
{
delete[] edge_data_1;
}
// edge 2
if (face.edge_orientation[2] == 1)
{
size_t idx = levelinfo::num_microvertices_per_face(level) + levelinfo::num_microedges_per_face(level) - 1;
for (size_t i = 0; i < rowsize; ++i)
{
// fmt::print("idx[{}] = {}\n", i, idx);
face_data[idx] = edge_data_2[i];
idx -= i+2;
}
idx = levelinfo::num_microvertices_per_face(level) + levelinfo::num_microedges_per_face(level) - 1 - 4;
for (size_t i = 1; i < rowsize-2; ++i)
{
// fmt::print("idx[{}] = {}\n", i, idx);
face_data[idx] = edge_data_2[edge_2_offset + i];
idx -= i+3;
}
}
else
{
size_t idx = 0;
for (size_t i = 0; i < rowsize; ++i)
{
// fmt::print("idx[{}] = {}\n", i, idx);
face_data[idx] = edge_data_2[i];
idx += rowsize-i;
}
idx = rowsize + 1;
for (size_t i = 1; i < rowsize-2; ++i)
{
// fmt::print("idx[{}] = {}\n", i, idx);
face_data[idx] = edge_data_2[edge_2_offset + i];
idx += rowsize-i;
}
}
if (face.edges[2]->rank != rk)
{
delete[] edge_data_2;
}
// else
// {
// size_t idx = levelinfo::num_microvertices_per_face(level) - 1;
// for (size_t i = 0; i < rowsize; ++i)
// {
// face_data[idx] = edge_data_1[i];
// idx -= i + 1;
// }
// }
// if (face.edges[1]->rank != rk)
// {
// delete[] edge_data_1;
// }
// // edge 2
// if (face.edge_orientation[2] == 1)
// {
// size_t idx = levelinfo::num_microvertices_per_face(level) - 1;
// for (size_t i = 0; i < rowsize; ++i)
// {
// face_data[idx] = edge_data_2[i];
// idx -= i+2;
// }
// }
// else
// {
// size_t idx = 0;
// for (size_t i = 0; i < rowsize; ++i)
// {
// face_data[idx] = edge_data_2[i];
// idx += rowsize-i;
// }
// }
// if (face.edges[2]->rank != rk)
// {
// delete[] edge_data_2;
// }
}
}
......
include/p2functionspace/p2function.hpp
View file @ 6e849737
......@@ -129,7 +129,7 @@ public:
{
if (edge.rank == rank)
{
P2Edge::interpolate(edge, memory_id, expr, level, testFlag(edge.type, flag));
P2Edge::interpolate(edge, memory_id, expr, level, flag);
}
}
......@@ -142,7 +142,7 @@ public:
{
if (face.rank == rank && testFlag(face.type, flag))
{
// P2Face::interpolate(face, memory_id, expr, level);
P2Face::interpolate(face, memory_id, expr, level);
}
}
}
......
test/test.cpp
View file @ 6e849737
......@@ -10,12 +10,12 @@ int main(int argc, char* argv[])
{
fmt::printf("TinyHHG default test\n");
}
hhg::Mesh mesh("../data/meshes/tri_1el.msh");
hhg::Mesh mesh("../data/meshes/bfs_126el.msh");
size_t minLevel = 2;
size_t maxLevel = 2;
size_t maxLevel = 4;
std::function<double(const hhg::Point3D&)> expr = [](const hhg::Point3D& x) { return 1.0; };
std::function<double(const hhg::Point3D&)> expr = [](const hhg::Point3D& x) { return cos(6*M_PI*x[0]); };
hhg::P2Function u("u", mesh, minLevel, maxLevel);
// hhg::P1Function Lu("Lu", mesh, minLevel, maxLevel);
......@@ -24,9 +24,9 @@ int main(int argc, char* argv[])
// hhg::P1MassOperator L_gen(mesh, minLevel, maxLevel);
hhg::VTKWriter({&u}, maxLevel, "../output", "test");
hhg::P2Edge::print(mesh.edges[4],u.memory_id,maxLevel);
hhg::P2Vertex::print(mesh.vertices[0], u.memory_id, maxLevel);
// hhg::P2Edge::print(mesh.edges[4],u.memory_id,maxLevel);
// hhg::P2Vertex::print(mesh.vertices[0], u.memory_id, maxLevel);
// hhg::P2Face::print(mesh.faces[1], u.memory_id, maxLevel);
MPI_Finalize();
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment