40 unsigned numElementsY,
41 unsigned numNodesPerEdge,
42 double proportionalGap,
47 assert(numElementsX > 0);
48 assert(numElementsY > 0);
49 assert(numNodesPerEdge > 2);
50 assert(proportionalGap > 0.0);
51 assert(padding > 0.0 && padding < 0.5);
54 unit_vector<double> x_unit(2,0);
55 unit_vector<double> y_unit(2,1);
57 std::vector<Node<2>*> nodes;
58 std::vector<ImmersedBoundaryElement<2,2>*> elements;
60 double width = 0.5 + 1.5 * numElementsX;
61 double height = sqrt(3.0) * numElementsY + 0.5 * sqrt(3.0) * (numElementsY > 1);
63 double max = width > height ? width : height;
65 double radius = (1.0 - 2.0 * padding) / max;
72 std::vector<c_vector<double, 2> > offsets(numElementsY * numElementsX);
75 c_vector<double, 2> global_offset;
78 global_offset[0] = padding + radius;
79 global_offset[1] = 0.5 - 0.5 * height + 0.5 * radius * sqrt(3.0);
83 global_offset[0] = 0.5 - 0.5 * width + radius;
84 global_offset[1] = padding + 0.5 * sqrt(3.0) * radius;
88 for (
unsigned x = 0; x < numElementsX; x++)
90 for (
unsigned y = 0; y < numElementsY; y++)
92 unsigned idx = x * numElementsY + y;
94 offsets[idx][0] = global_offset[0] + 1.5 * radius * x;
95 offsets[idx][1] = global_offset[1] + 0.5 * sqrt(3.0) * radius * (2.0 * y + (
double)(x % 2 == 1));
100 std::vector<c_vector<double, 2> > node_locations =
GetUnitHexagon(numNodesPerEdge);
102 double scale = 1.0 - proportionalGap;
105 for (
unsigned offset = 0; offset < offsets.size(); offset++)
107 std::vector<Node<2>*> nodes_this_elem;
109 for (
unsigned location = 0; location < node_locations.size(); location++)
111 unsigned index = offset * node_locations.size() + location;
112 Node<2>* p_node =
new Node<2>(index, offsets[offset] + scale * radius * node_locations[location],
true);
114 nodes_this_elem.push_back(p_node);
115 nodes.push_back(p_node);
121 if (offset < numElementsY ||
122 offset >= numElementsY * (numElementsX - 1) ||
123 offset % numElementsY == 0 ||
124 (offset + 1) % numElementsY == 0)
129 elements.push_back(p_elem);
148 std::vector<c_vector<double, 2> > locations(numPtsPerSide * 6);
151 std::vector<c_vector<double, 2> > vertices(7);
152 double sixty_degrees = M_PI / 3.0;
153 for (
unsigned vertex = 0; vertex < vertices.size(); vertex++)
155 vertices[vertex][0] = cos((
double)vertex * sixty_degrees);
156 vertices[vertex][1] = sin((
double)vertex * sixty_degrees);
160 for (
unsigned vertex = 0; vertex < 6; vertex++)
162 c_vector<double, 2> this_vertex = vertices[vertex];
163 c_vector<double, 2> next_vertex = vertices[vertex + 1];
164 c_vector<double, 2> vec_between = next_vertex - this_vertex;
166 for (
unsigned i = 0; i < numPtsPerSide; i++)
168 locations[vertex * numPtsPerSide + i] = this_vertex + i * vec_between / (
double)numPtsPerSide;