Commit 8879de3d authored by abbasaa's avatar abbasaa
Browse files

refactor and comments

parent efd9bde4
......@@ -6,10 +6,9 @@
#include <fstream>
#include <iomanip>
#include <limits>
#include <getopt.h>
#include <stdio.h>
#include "face.h"
#include "face.h"
using namespace std;
......@@ -114,8 +113,6 @@ int boxGenMain(map<string,string>& config, string protein_filename) {
getProteinLimits(protein_filename, &protein_min_x, &protein_max_x, &protein_min_y, &protein_max_y, &membrane_bottom_height, &membrane_top_height);
membrane_bottom_height = stod(getVal(config, "membrane_bottom_height"));
membrane_top_height = stod(getVal(config, "membrane_top_height"));
cout << "membrane_bottom_height = " << to_string(membrane_bottom_height) << endl;
cout << "membrane_top_height = " << to_string(membrane_top_height) << endl;
int nz_middle = int(membrane_top_height - membrane_bottom_height);
......@@ -265,7 +262,7 @@ int boxGenMain(map<string,string>& config, string protein_filename) {
}
n_nodes += add_nodes.size();
// FILE OUTPUT
// File output
poly_file << "#POLY FILE OUTPUT\n";
poly_file << (n_nodes + protein_nodes) << "\t3\t0\t0\n";
poly_file << node_out.str();
......
......@@ -91,26 +91,26 @@ public:
string node_file_name = "";
string ele_file_name = "";
string neigh_file_name = "";
string outFileName = "";
string out_file_name = "";
// Minimum and maximum z coordinates of the membrane
double z_min = -11, z_max = 6;
double z_min = 0, z_max = 0;
private:
// Store properties and component of the volume mesh
double x_max = -std::numeric_limits<double>::infinity();
double y_max = -std::numeric_limits<double>::infinity();
int num_nodes = 0, num_tetras = 0;
vector<Mem_Node> nodes;
vector<Tetra> tetrahedrons;
vector<Neigh> neighbors;
vector<int> region_markers;
vector<bool> visited;
int numNodes = 0, numTetras = 0;
double x_max = -std::numeric_limits<double>::infinity();
double y_max = -std::numeric_limits<double>::infinity();
// Vectors and values used in membrane extraction BFS
int tetra_index = -1;
vector<int> region_markers;
vector<bool> visited;
int protein_marker = 2, solvent_marker = 1, mem_marker = 3;
// Markers used for regions in output.xml file
int solvent_marker = 1, protein_marker = 2, mem_marker = 3;
};
......@@ -132,7 +132,7 @@ int memGenMain(map<string,string>& config, string file_prefix) {
generator.node_file_name = file_prefix + ".node";
generator.ele_file_name = file_prefix + ".ele";
generator.neigh_file_name = file_prefix + ".neigh";
generator.outFileName = "output.xml";
generator.out_file_name = "output.xml";
generator.z_min = stod(getVal(config, "membrane_bottom_height"));
generator.z_max = stod(getVal(config, "membrane_top_height"));
......@@ -156,17 +156,17 @@ void MemGen::readFiles() {
int temp;
// Read in number of nodes and tetrahedra in mesh
node_file >> numNodes >> temp >> temp >> temp;
ele_file >> numTetras >> temp >> temp;
node_file >> num_nodes >> temp >> temp >> temp;
ele_file >> num_tetras >> temp >> temp;
nodes.reserve(numNodes);
tetrahedrons.reserve(numTetras);
neighbors.reserve(numTetras);
region_markers.resize(numTetras);
visited.resize(numTetras, false);
nodes.reserve(num_nodes);
tetrahedrons.reserve(num_tetras);
neighbors.reserve(num_tetras);
region_markers.resize(num_tetras);
visited.resize(num_tetras, false);
cout << "Reading .node file\n";
for (int i = 0; i < numNodes; i++){
for (int i = 0; i < num_nodes; i++){
double x, y, z;
node_file >> temp >> x >> y >> z;
x_max = (x > x_max) ? x : x_max;
......@@ -176,7 +176,7 @@ void MemGen::readFiles() {
cout << "Reading .ele file\n";
set<int> markers;
for (int i = 0; i < numTetras; ++i){
for (int i = 0; i < num_tetras; ++i){
int node1, node2, node3, node4, marker;
ele_file >> temp >> node1 >> node2 >> node3 >> node4 >> marker;
node1 -= 1;
......@@ -196,7 +196,7 @@ void MemGen::readFiles() {
cout << "Reading .neigh file\n";
neigh_file >> temp >> temp;
for (int i = 0; i < numTetras; ++i){
for (int i = 0; i < num_tetras; ++i){
int neigh1, neigh2, neigh3, neigh4;
neigh_file >> temp >> neigh1 >> neigh2 >> neigh3 >> neigh4;
neigh1 -= 1; // zero-indexes the tetrahedra in the .ele file
......@@ -217,7 +217,7 @@ void MemGen::readFiles() {
void MemGen::findStartingTetra() {
cout << "Finding point's tetrahedron. NumTetras: " << numTetras << "\n";
cout << "Finding point's tetrahedron. num_tetras: " << num_tetras << "\n";
double x_mem = x_max - 1.11;
double y_mem = y_max - 1.11;
double z_mem = z_max - 1.11;
......@@ -225,7 +225,7 @@ void MemGen::findStartingTetra() {
// Iterate through all tetrahedra in mesh until
// tetrahedra containing point p is found
for (int i = 0; i < numTetras; ++i){
for (int i = 0; i < num_tetras; ++i){
const Mem_Node& v1 = nodes[tetrahedrons[i].node1];
const Mem_Node& v2 = nodes[tetrahedrons[i].node2];
const Mem_Node& v3 = nodes[tetrahedrons[i].node3];
......@@ -272,29 +272,29 @@ void MemGen::extractMembrane() {
void MemGen::genOutput() {
ofstream outFile;
outFile.open(outFileName);
ofstream out_file;
out_file.open(out_file_name);
// File Output
cout << "Outputting to file " << outFileName << "\n";
outFile << "<?xml version=\"1.0\"?>\n";
outFile << "<dolfin xmlns:dolfin=\"http://fenicsproject.org\">\n";
outFile << "\t<mesh celltype=\"tetrahedron\" dim=\"3\">\n";
outFile << "\t\t<vertices size=\"" << numNodes << "\">\n";
for (int i = 0; i < numNodes; i++){
outFile << "\t\t\t<vertex index=\"" << i << "\" x=\"" << nodes[i].x << "\" y=\""
cout << "Outputting to file " << out_file_name << "\n";
out_file << "<?xml version=\"1.0\"?>\n";
out_file << "<dolfin xmlns:dolfin=\"http://fenicsproject.org\">\n";
out_file << "\t<mesh celltype=\"tetrahedron\" dim=\"3\">\n";
out_file << "\t\t<vertices size=\"" << num_nodes << "\">\n";
for (int i = 0; i < num_nodes; i++){
out_file << "\t\t\t<vertex index=\"" << i << "\" x=\"" << nodes[i].x << "\" y=\""
<< nodes[i].y << "\" z=\"" << nodes[i].z << "\" />\n";
}
outFile << "\t\t</vertices>\n";
outFile << "\t\t<cells size=\"" << numTetras << "\">\n";
for(int i = 0; i < numTetras; ++i){
outFile << "\t\t\t<tetrahedron index=\"" << i << "\" v0=\"" << tetrahedrons[i].node1
out_file << "\t\t</vertices>\n";
out_file << "\t\t<cells size=\"" << num_tetras << "\">\n";
for(int i = 0; i < num_tetras; ++i){
out_file << "\t\t\t<tetrahedron index=\"" << i << "\" v0=\"" << tetrahedrons[i].node1
<< "\" v1=\"" << tetrahedrons[i].node2 << "\" v2=\"" << tetrahedrons[i].node3
<< "\" v3=\"" << tetrahedrons[i].node4 << "\" />\n";
}
outFile << "\t\t</cells>\n\t\t<data>\n\t\t\t<data_entry name=\"cell domains\">\n";
outFile << "\t\t\t\t<mesh_function type=\"uint\" dim=\"3\" size=\"" << numTetras << "\">\n";
for (int i = 0; i < numTetras; ++i){
out_file << "\t\t</cells>\n\t\t<data>\n\t\t\t<data_entry name=\"cell domains\">\n";
out_file << "\t\t\t\t<mesh_function type=\"uint\" dim=\"3\" size=\"" << num_tetras << "\">\n";
for (int i = 0; i < num_tetras; ++i){
int out_marker;
if (region_markers[i] == solvent_marker)
out_marker = 1;
......@@ -304,11 +304,11 @@ void MemGen::genOutput() {
out_marker = 2;
}
outFile << "\t\t\t\t\t<entity index=\"" << i << "\" value=\"" << out_marker << "\" />\n";
out_file << "\t\t\t\t\t<entity index=\"" << i << "\" value=\"" << out_marker << "\" />\n";
}
outFile << "\t\t\t\t</mesh_function>\n\t\t\t</data_entry>\n\t\t</data>\n\t</mesh>\n</dolfin>" << endl;
out_file << "\t\t\t\t</mesh_function>\n\t\t\t</data_entry>\n\t\t</data>\n\t</mesh>\n</dolfin>" << endl;
outFile.close();
out_file.close();
}
bool MemGen::inMembraneZ(const Tetra& tet){
......
from fenics import *
#from pathlib import Path
import os
import glob
def main():
mesh = Mesh('output.xml') # Get the mesh of a box domain
mesh = Mesh('output.xml') # Read generated mesh
# cell_domains: an object of MeshFunction defined on the set of tetrahedrons
# since mesh.topology().dim() is 3
# since mesh.topology().dim() is 3
cell_domains = MeshFunction('size_t', mesh, mesh.topology().dim())
# Get the labels of tetrahedrons from the mesh data with the key word "cell domains"
cell_domains.array()[:] = mesh.data().array("cell domains", mesh.topology().dim())
DsMesh = SubMesh(mesh, cell_domains, 1) # Mesh of the solvent region Ds
DpMesh = SubMesh(mesh, cell_domains, 2) # Mesh of the protein region Dp
DmMesh = SubMesh(mesh, cell_domains, 3) # Mesh of the solvent region Ds
DmMesh = SubMesh(mesh, cell_domains, 3) # Mesh of the membrane region Dm
# Save submesh in .pvd format
File('protein.pvd') << DpMesh
File('solvent.pvd') << DsMesh
File('membrane.pvd') << DmMesh
# Then you can display them by paraview
File('../output/protein.pvd') << DpMesh
File('../output/solvent.pvd') << DsMesh
File('../output/membrane.pvd') << DmMesh
# Remove extra files created by program
for TMSmesh_file in glob.glob('pdb*'):
print("Deleting {}".format(TMSmesh_file))
os.remove(TMSmesh_file)
other_files = [ "output.xml" ]
for del_file in other_files:
print("Deleting {}".format(del_file))
os.remove(del_file)
if __name__ == "__main__":
main()
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