37 #include "HoneycombMeshGenerator.hpp" 39 #include <boost/foreach.hpp> 40 #include "TrianglesMeshReader.hpp" 41 #include "OutputFileHandler.hpp" 42 #include "RandomNumberGenerator.hpp" 48 mMeshFilename(
"mesh"),
49 mDomainWidth(numNodesAlongWidth*scaleFactor),
50 mNumCellWidth(numNodesAlongWidth),
51 mNumCellLength(numNodesAlongLength)
57 assert(scaleFactor > 0.0);
60 std::stringstream pid;
62 OutputFileHandler output_file_handler(
"2D_temporary_honeycomb_mesh_" + pid.str());
68 double vertical_spacing = (sqrt(3.0)/2)*horizontal_spacing;
74 num_nodes_along_width += 2*ghosts;
75 num_nodes_along_length += 2*ghosts;
77 unsigned num_nodes = num_nodes_along_width*num_nodes_along_length;
78 unsigned num_elem_along_width = num_nodes_along_width-1;
79 unsigned num_elem_along_length = num_nodes_along_length-1;
80 unsigned num_elem = 2*num_elem_along_width*num_elem_along_length;
81 unsigned num_edges = 3*num_elem_along_width*num_elem_along_length + num_elem_along_width + num_elem_along_length;
83 double x0 = -horizontal_spacing*ghosts;
84 double y0 = -vertical_spacing*ghosts;
86 mBottom = -vertical_spacing*ghosts;
87 mTop =
mBottom + vertical_spacing*(num_nodes_along_length-1);
90 out_stream p_node_file = output_file_handler.OpenOutputFile(
mMeshFilename+
".node");
91 (*p_node_file) << std::scientific;
92 (*p_node_file) << std::setprecision(20);
93 (*p_node_file) << num_nodes <<
"\t2\t0\t1" << std::endl;
96 for (
unsigned i=0; i<num_nodes_along_length; i++)
98 for (
unsigned j=0; j<num_nodes_along_width; j++)
108 unsigned boundary = 0;
109 if (i==0 || i==num_nodes_along_length-1 || j==0 || j==num_nodes_along_width-1)
114 double x = x0 + horizontal_spacing*((
double)j + 0.25*(1.0+
SmallPow(-1.0,i+1)));
115 double y = y0 + vertical_spacing*(
double)i;
118 if ((y<0.0) && (y>-1e-12))
126 (*p_node_file) << node++ <<
"\t" << x <<
"\t" << y <<
"\t" << boundary << std::endl;
129 p_node_file->close();
132 out_stream p_elem_file = output_file_handler.OpenOutputFile(
mMeshFilename+
".ele");
133 (*p_elem_file) << std::scientific;
135 out_stream p_edge_file = output_file_handler.OpenOutputFile(
mMeshFilename+
".edge");
136 (*p_node_file) << std::scientific;
138 (*p_elem_file) << num_elem <<
"\t3\t0" << std::endl;
139 (*p_edge_file) << num_edges <<
"\t1" << std::endl;
143 for (
unsigned i=0; i<num_elem_along_length; i++)
145 for (
unsigned j=0; j < num_elem_along_width; j++)
147 unsigned node0 = i*num_nodes_along_width + j;
148 unsigned node1 = i*num_nodes_along_width + j+1;
149 unsigned node2 = (i+1)*num_nodes_along_width + j;
156 (*p_elem_file) << elem++ <<
"\t" << node0 <<
"\t" << node1 <<
"\t" << node2 << std::endl;
158 unsigned horizontal_edge_is_boundary_edge = 0;
159 unsigned vertical_edge_is_boundary_edge = 0;
162 horizontal_edge_is_boundary_edge = 1;
166 vertical_edge_is_boundary_edge = 1;
169 (*p_edge_file) << edge++ <<
"\t" << node0 <<
"\t" << node1 <<
"\t" << horizontal_edge_is_boundary_edge << std::endl;
170 (*p_edge_file) << edge++ <<
"\t" << node1 <<
"\t" << node2 <<
"\t" << 0 << std::endl;
171 (*p_edge_file) << edge++ <<
"\t" << node2 <<
"\t" << node0 <<
"\t" << vertical_edge_is_boundary_edge << std::endl;
173 node0 = i*num_nodes_along_width + j + 1;
179 node1 = (i+1)*num_nodes_along_width + j+1;
180 node2 = (i+1)*num_nodes_along_width + j;
182 (*p_elem_file) << elem++ <<
"\t" << node0 <<
"\t" << node1 <<
"\t" << node2 << std::endl;
186 for (
unsigned i=0; i<num_elem_along_length; i++)
188 unsigned node0, node1;
192 node0 = (i+1)*num_nodes_along_width - 1;
193 node1 = (i+2)*num_nodes_along_width - 1;
197 node0 = (i+1)*num_nodes_along_width;
198 node1 = (i)*num_nodes_along_width;
200 (*p_edge_file) << edge++ <<
"\t" << node0 <<
"\t" << node1 <<
"\t" << 1 << std::endl;
203 for (
unsigned j=0; j<num_elem_along_width; j++)
205 unsigned node0 = num_nodes_along_width*(num_nodes_along_length-1) + j;
206 unsigned node1 = num_nodes_along_width*(num_nodes_along_length-1) + j+1;
207 (*p_edge_file) << edge++ <<
"\t" << node1 <<
"\t" << node0 <<
"\t" << 1 << std::endl;
210 p_elem_file->close();
211 p_edge_file->close();
222 output_file_handler.FindFile(
"").Remove();
240 std::vector<unsigned> location_indices;
246 location_indices.push_back(i);
249 return location_indices;
256 EXCEPTION(
"Cannot call GetCircularMesh on a HoneycombMesh with ghost nodes");
260 c_vector<double,2> centre = zero_vector<double>(2);
279 c_vector<double,2>& r_location =
mpMesh->
GetNode(i)->rGetModifiableLocation();
280 c_vector<double,2> shift;
282 double max_jiggle = radius*5e-6;
283 shift[0] = max_jiggle*(p_gen->
ranf()-0.5);
284 shift[1] = max_jiggle*(p_gen->
ranf()-0.5);
std::string mMeshFilename
double SmallPow(double x, unsigned exponent)
void DeleteNodePriorToReMesh(unsigned index)
MutableMesh< 2, 2 > * GetCircularMesh(double radius)
virtual void Translate(const c_vector< double, SPACE_DIM > &rDisplacement)
Node< SPACE_DIM > * GetNode(unsigned index) const
#define EXCEPTION(message)
void ConstructFromMeshReader(AbstractMeshReader< ELEMENT_DIM, SPACE_DIM > &rMeshReader)
std::string GetOutputDirectoryFullPath() const
virtual unsigned GetNumAllNodes() const
unsigned GetNumNodes() const
std::vector< unsigned > GetCellLocationIndices()
std::set< unsigned > mGhostNodeIndices
virtual void ReMesh(NodeMap &map)
static RandomNumberGenerator * Instance()
virtual MutableMesh< 2, 2 > * GetMesh()
void SetMeshHasChangedSinceLoading()
virtual ~HoneycombMeshGenerator()
MutableMesh< 2, 2 > * mpMesh