43 std::vector<CellPtr>& rCells,
44 const std::vector<unsigned> locationIndices,
46 double ghostCellSpringStiffness,
47 double ghostGhostSpringStiffness,
48 double ghostSpringRestLength)
50 mGhostCellSpringStiffness(ghostCellSpringStiffness),
51 mGhostGhostSpringStiffness(ghostGhostSpringStiffness),
52 mGhostSpringRestLength(ghostSpringRestLength)
54 if (!locationIndices.empty())
57 std::set<unsigned> node_indices;
58 std::set<unsigned> location_indices;
59 std::set<unsigned> ghost_node_indices;
63 node_indices.insert(this->
GetNode(i)->GetIndex());
65 for (
unsigned i=0; i<locationIndices.size(); i++)
67 location_indices.insert(locationIndices[i]);
70 std::set_difference(node_indices.begin(), node_indices.end(),
71 location_indices.begin(), location_indices.end(),
72 std::inserter(ghost_node_indices, ghost_node_indices.begin()));
152 assert(rNodeAGlobalIndex != rNodeBGlobalIndex);
153 c_vector<double, DIM> unit_difference;
154 const c_vector<double, DIM>& r_node_a_location = this->GetNode(rNodeAGlobalIndex)->rGetLocation();
155 const c_vector<double, DIM>& r_node_b_location = this->GetNode(rNodeBGlobalIndex)->rGetLocation();
158 unit_difference = this->mrMesh.GetVectorFromAtoB(r_node_a_location, r_node_b_location);
160 double distance_between_nodes = norm_2(unit_difference);
161 unit_difference /= distance_between_nodes;
163 double rest_length = 1.0;
164 double spring_stiffness = mGhostCellSpringStiffness;
165 if (this->mIsGhostNode[rNodeAGlobalIndex] && this->mIsGhostNode[rNodeBGlobalIndex])
167 rest_length = mGhostSpringRestLength;
168 spring_stiffness = mGhostGhostSpringStiffness;
171 return spring_stiffness * unit_difference * (distance_between_nodes - rest_length);
300 std::vector<c_vector<double, DIM> > drdt(this->GetNumNodes());
301 for (
unsigned i=0; i<drdt.size(); i++)
303 drdt[i] = zero_vector<double>(DIM);
311 unsigned nodeA_global_index = edge_iterator.GetNodeA()->GetIndex();
312 unsigned nodeB_global_index = edge_iterator.GetNodeB()->GetIndex();
314 c_vector<double, DIM> force = CalculateForceBetweenGhostNodes(nodeA_global_index, nodeB_global_index);
316 if (!this->mIsGhostNode[nodeA_global_index])
318 drdt[nodeB_global_index] -= force;
322 drdt[nodeA_global_index] += force;
324 if (this->mIsGhostNode[nodeB_global_index])
326 drdt[nodeB_global_index] -= force;
332 node_iter != this->mrMesh.GetNodeIteratorEnd();
335 unsigned node_index = node_iter->GetIndex();
336 if (this->mIsGhostNode[node_index])
338 node_iter->ClearAppliedForce();
339 node_iter->AddAppliedForceContribution(drdt[node_index]);
364 std::stringstream time;
365 time << num_timesteps;
367 if (this->mWriteVtkAsPoints)
373 unsigned num_vtk_cells = this->rGetMesh().GetNumNodes();
375 cell_writer_iter != this->mCellWriters.end();
379 std::vector<double> vtk_cell_data(num_vtk_cells);
383 node_iter != this->rGetMesh().GetNodeIteratorEnd();
387 unsigned node_index = node_iter->GetIndex();
390 if (this->IsGhostNode(node_index))
393 vtk_cell_data[node_index] = -1.0;
398 CellPtr p_cell = this->GetCellUsingLocationIndex(node_index);
401 vtk_cell_data[node_index] = (*cell_writer_iter)->GetCellDataForVtkOutput(p_cell,
this);
405 mesh_writer.
AddPointData((*cell_writer_iter)->GetVtkCellDataName(), vtk_cell_data);
410 std::vector<double> ghosts(num_vtk_cells);
412 node_iter != this->rGetMesh().GetNodeIteratorEnd();
415 unsigned node_index = node_iter->GetIndex();
416 ghosts[node_index] = (
double) (this->IsGhostNode(node_index));
423 *(this->mpVtkMetaFile) <<
" <DataSet timestep=\"";
424 *(this->mpVtkMetaFile) << num_timesteps;
425 *(this->mpVtkMetaFile) <<
"\" group=\"\" part=\"0\" file=\"mesh_results_";
426 *(this->mpVtkMetaFile) << num_timesteps;
427 *(this->mpVtkMetaFile) <<
".vtu\"/>\n";
429 if (this->mpVoronoiTessellation !=
nullptr)
435 unsigned num_vtk_cells = this->mpVoronoiTessellation->GetNumElements();
437 cell_writer_iter != this->mCellWriters.end();
441 std::vector<double> vtk_cell_data(num_vtk_cells);
445 elem_iter != this->mpVoronoiTessellation->GetElementIteratorEnd();
449 unsigned elem_index = elem_iter->GetIndex();
450 unsigned node_index = this->mpVoronoiTessellation->GetDelaunayNodeIndexCorrespondingToVoronoiElementIndex(elem_index);
453 if (this->IsGhostNode(node_index))
456 vtk_cell_data[elem_index] = -1.0;
461 CellPtr p_cell = this->GetCellUsingLocationIndex(node_index);
464 vtk_cell_data[elem_index] = (*cell_writer_iter)->GetCellDataForVtkOutput(p_cell,
this);
468 mesh_writer.
AddCellData((*cell_writer_iter)->GetVtkCellDataName(), vtk_cell_data);
473 std::vector<double> ghosts(num_vtk_cells);
475 elem_iter != this->mpVoronoiTessellation->GetElementIteratorEnd();
478 unsigned elem_index = elem_iter->GetIndex();
479 unsigned node_index = this->mpVoronoiTessellation->GetDelaunayNodeIndexCorrespondingToVoronoiElementIndex(elem_index);
480 ghosts[elem_index] = (
double) (this->IsGhostNode(node_index));
487 *(this->mpVtkMetaFile) <<
" <DataSet timestep=\"";
488 *(this->mpVtkMetaFile) << num_timesteps;
489 *(this->mpVtkMetaFile) <<
"\" group=\"\" part=\"0\" file=\"voronoi_results_";
490 *(this->mpVtkMetaFile) << num_timesteps;
491 *(this->mpVtkMetaFile) <<
".vtu\"/>\n";