Chaste Commit::ca8ccdedf819b6e02855bc0e8e6f50bdecbc5208
HoneycombVertexMeshGenerator.cpp
1/*
2
3Copyright (c) 2005-2024, University of Oxford.
4All rights reserved.
5
6University of Oxford means the Chancellor, Masters and Scholars of the
7University of Oxford, having an administrative office at Wellington
8Square, Oxford OX1 2JD, UK.
9
10This file is part of Chaste.
11
12Redistribution and use in source and binary forms, with or without
13modification, are permitted provided that the following conditions are met:
14 * Redistributions of source code must retain the above copyright notice,
15 this list of conditions and the following disclaimer.
16 * Redistributions in binary form must reproduce the above copyright notice,
17 this list of conditions and the following disclaimer in the documentation
18 and/or other materials provided with the distribution.
19 * Neither the name of the University of Oxford nor the names of its
20 contributors may be used to endorse or promote products derived from this
21 software without specific prior written permission.
22
23THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33
34*/
35
36#include "HoneycombVertexMeshGenerator.hpp"
37#include <boost/make_shared.hpp>
38#include <boost/shared_ptr.hpp>
39
41 unsigned numElementsUp,
42 bool isFlatBottom,
43 double cellRearrangementThreshold,
44 double t2Threshold,
45 double elementArea)
46{
47 assert(numElementsAcross > 0);
48 assert(numElementsUp > 0);
49 assert(cellRearrangementThreshold > 0.0);
50 assert(t2Threshold > 0.0);
51 assert(elementArea > 0.0);
52
53 std::vector<Node<2>*> nodes;
54 std::vector<VertexElement<2,2>*> elements;
55
56 unsigned node_index = 0;
57 unsigned node_indices[6];
58 unsigned element_index;
59
60 /*
61 * Create the nodes, row by row, from the bottom up. On the first row we
62 * have numElementsAcross nodes, all of which are boundary nodes. This row
63 * only exists if isFlatBottom is false.
64 */
65 if (!isFlatBottom)
66 {
67 for (unsigned i = 0; i < numElementsAcross; i++)
68 {
69 Node<2>* p_node = new Node<2>(node_index, true, i + 0.5, 0);
70 nodes.push_back(p_node);
71 node_index++;
72 }
73 }
74 /*
75 * On each interior row, and the first row if isFlatBottom is true, we have
76 * numElementsAcross+1 nodes. On the penultimate row, all nodes are boundary
77 * nodes. The same is true for the second row, if isFlatBottom is false. On
78 * other rows, the first and last nodes only are boundary nodes.
79 */
80 for (unsigned j = 1; j < 2*numElementsUp + 1; j++)
81 {
82 if ((j == 1) && isFlatBottom)
83 {
84 for (unsigned i = 0; i <= numElementsAcross; i++)
85 {
86 double x_coord = ((j%4 == 0)||(j%4 == 3)) ? i + 0.5 : i;
87 double y_coord = 0.0;
88
89 Node<2>* p_node = new Node<2>(node_index, true, x_coord, y_coord);
90 nodes.push_back(p_node);
91 node_index++;
92 }
93 }
94 else
95 {
96 for (unsigned i = 0; i <= numElementsAcross; i++)
97 {
98 double x_coord = ((j%4 == 0)||(j%4 == 3)) ? i + 0.5 : i;
99 double y_coord = (1.5*j - 0.5*(j%2))*0.5/sqrt(3.0);
100 bool is_boundary_node = (j==1 || j==2*numElementsUp || i==0 || i==numElementsAcross) ? true : false;
101
102 Node<2>* p_node = new Node<2>(node_index, is_boundary_node, x_coord, y_coord);
103 nodes.push_back(p_node);
104 node_index++;
105 }
106 }
107 }
108
109 /*
110 * On the last row we have numElementsAcross nodes, all of which are boundary nodes.
111 */
112 double y_coord = (1.5*(2*numElementsUp+1) - 0.5*((2*numElementsUp+1)%2))*0.5/sqrt(3.0);
113 if (((2*numElementsUp+1)%4 == 0)||((2*numElementsUp+1)%4 == 3))
114 {
115 Node<2>* p_node = new Node<2>(node_index, true, 0.5, y_coord);
116 nodes.push_back(p_node);
117 node_index++;
118 }
119 for (unsigned i=1; i<numElementsAcross; i++)
120 {
121 double x_coord = (((2*numElementsUp+1)%4 == 0)||((2*numElementsUp+1)%4 == 3)) ? i+0.5 : i;
122
123 Node<2>* p_node = new Node<2>(node_index, true, x_coord, y_coord);
124 nodes.push_back(p_node);
125 node_index++;
126 }
127 if (((2*numElementsUp+1)%4 == 1)||((2*numElementsUp+1)%4 == 2))
128 {
129 Node<2>* p_node = new Node<2>(node_index, true, numElementsAcross, y_coord);
130 nodes.push_back(p_node);
131 node_index++;
132 }
133
134 /*
135 * Create the elements. The array node_indices contains the
136 * global node indices from bottom, going anticlockwise.
137 */
138 for (unsigned j = 0; j < numElementsUp; j++)
139 {
140 if (j == 0 && isFlatBottom)
141 {
142 unsigned short_node_indices[5];
143 for (unsigned i = 0; i < numElementsAcross; i++)
144 {
145 short_node_indices[0] = i;
146 short_node_indices[1] = i + 1;
147 short_node_indices[2] = short_node_indices[1] + numElementsAcross + 1;
148 short_node_indices[3] = short_node_indices[2] + numElementsAcross;
149 short_node_indices[4] = short_node_indices[2] - 1;
150
151 std::vector<Node<2>*> element_nodes;
152 for (unsigned k = 0; k < 5; k++)
153 {
154 element_nodes.push_back(nodes[short_node_indices[k]]);
155 }
156
157 element_index = j*numElementsAcross + i;
158 VertexElement<2,2>* p_element = new VertexElement<2,2>(element_index, element_nodes);
159 elements.push_back(p_element);
160 }
161 }
162 else
163 {
164 for (unsigned i = 0; i < numElementsAcross; i++)
165 {
166 if (j == 0)
167 {
168 node_indices[0] = i;
169 }
170 else
171 {
172 unsigned first_index;
173 if (isFlatBottom)
174 {
175 first_index = 2*j*(numElementsAcross+1) - 1*(j%2==0) + i - numElementsAcross;
176 }
177 else
178 {
179 first_index = 2*j*(numElementsAcross+1) - 1*(j%2==0) + i;
180 }
181 node_indices[0] = first_index; // different for even/odd rows
182 }
183 node_indices[1] = node_indices[0] + numElementsAcross + 1 + 1*(j%2==0 && j>0);
184 node_indices[2] = node_indices[1] + numElementsAcross + 1;
185 node_indices[3] = node_indices[2] + numElementsAcross + 1*(j%2==1 && j<numElementsUp-1);
186 node_indices[4] = node_indices[2] - 1;
187 node_indices[5] = node_indices[1] - 1;
188
189 std::vector<Node<2>*> element_nodes;
190 for (unsigned k = 0; k < 6; k++)
191 {
192 element_nodes.push_back(nodes[node_indices[k]]);
193 }
194
195 element_index = j*numElementsAcross + i;
196 VertexElement<2,2>* p_element = new VertexElement<2,2>(element_index, element_nodes);
197 elements.push_back(p_element);
198 }
199 }
200 }
201
202 mpMesh = boost::make_shared<MutableVertexMesh<2,2> >(nodes, elements, cellRearrangementThreshold, t2Threshold);
203
204
205 // Scale the mesh so that each element's area takes the value elementArea
206 mpMesh->Scale(sqrt(elementArea*2.0/sqrt(3.0)), sqrt(elementArea*2.0/sqrt(3.0)));
207}
208
209boost::shared_ptr<MutableVertexMesh<2,2> > HoneycombVertexMeshGenerator::GetMesh()
210{
211 return mpMesh;
212}
boost::shared_ptr< MutableVertexMesh< 2, 2 > > mpMesh
virtual boost::shared_ptr< MutableVertexMesh< 2, 2 > > GetMesh()
Definition Node.hpp:59