CSMP code snippets
Here we will post various bits of CSMP++ code that you might require in the course of your studies. Clicking one the links below will reveal the actual code snippet.
bool El_vs_Group ( SuperGroup<csp_float,DIM>& sg, stl_index element_id, const char* group ) { bool out(false); std::map<stl_index,DATA_STYLE>& GrEl = sg.ReferenceGroup(group).ElementIDMap(); std::map<stl_index,DATA_STYLE>::iterator it = GrEl.find(element_id); if ( it != GrEl.end() ) out = true; return out; }
for ( std::map<stl_index,DATA_STYLE>::const_iterator it = sg.ReferenceGroup("GRAINS").NodesBegin(); it != sg.ReferenceGroup("GRAINS").NodesEnd(); it++ ) { if ( (*it).second == BOUNDARY ) { //your code } }
csp_float NodeVolume ( SuperGroup<csp_float,DIM>& sg, stl_index node_id ) { csp_float total_area = 0.; const MemoryManager<csp_float,DIM>& pmem(sg.ReferencePropertyStorage()); stl_index parents( sg.ReferenceMesh().N(node_id-1U).Parents() ); for( stl_index parent = 0U; parent < parents; parent++ ) { const stl_index parent_id(sg.ReferenceMesh().N(node_id-1U).ParentElementID(parent)); const Element<csp_float,DIM>& elmt(sg.ReferenceMesh().E(parent_id-1U)); ElementFiniteVolumeTraits<csp_float,DIM,Element<csp_float,DIM>,MemoryManager<csp_float,DIM> > efvt(elmt,pmem); total_area += efvt.SectorVolume(sg.ReferenceMesh().N(node_id-1U).ParentLocalNodeNumber(parent)); } return total_area; }
// Create a face variable to turn faces on. // Then use this command to actually set the status of boundary faces to BOUNDARY sg.ChangePropertyStatusInGroupToWhere( "GRAINS", "face var", BOUNDARY, BOUNDARY ); csp_float BoundaryNodeSurface ( SuperGroup<csp_float,2>& sg, stl_index node_id, const char* variable ) { csp_float total_area = 0.0; std::vector<stl_index> nd_ids, fc_ids; bool hasNode, isBoundary, isSGBoundary; const MemoryManager<csp_float,2>& pmem(sg.ReferencePropertyStorage()); stl_index parents( sg.ReferenceMesh().N(node_id-1U).Parents() ); for( stl_index parent = 0U; parent < parents; parent++ ) { stl_index parent_id(sg.ReferenceMesh().N(node_id-1U).ParentElementID(parent)); const Element<csp_float,2>& elmt(sg.ReferenceMesh().E(parent_id-1U)); elmt.FaceIDs( fc_ids ); stl_index faces( elmt.Faces() ); for( stl_index face = 0U; face < faces; face++ ) { const Face<csp_float,2>& fc(sg.ReferenceMesh().F(fc_ids[face]-1U)); fc.NodeIDs( nd_ids ); hasNode = (std::find(nd_ids.begin(), nd_ids.end(), node_id) != nd_ids.end()); isBoundary = (fc.Status(pmem, sg.ReferencePropertyDatabase().StorageKey(variable)) == BOUNDARY); isSGBoundary = (fc.AtBoundary() != NOT); if ( hasNode and isBoundary ) { if ( isSGBoundary ) total_area += fc.Area()/2.0; else total_area += fc.Area()/4.0; } } } return total_area; }
void scaleRegion( SuperGroup<csp_float,2>& sg, csp_float scale_factor) { csp_float x_, y_; for ( std::deque<Node<csp_float,2> >::iterator it = sg.ReferenceMesh().NodesBegin(); it != sg.ReferenceMesh().NodesEnd(); it++ ) { x_ = it->x(); it->x(x_/scale_factor); y_ = it->y(); it->y(y_/scale_factor); } }
// The curves are at 4 equidistant positions of the model, template<typename fT, stl_index dim> void monitorBTC( SuperGroup<fT,dim>& sg, fT flux ) { static bool first_time(true); static fT x[4]; static fT vol[4]; static fT time0; static typename std::vector<std::list<stl_index> > eid(4); static typename std::vector<std::list<fT> > vel(4); static typename std::vector<std::ofstream*> ofs(4); static csp::Index c_key(sg.ReferencePropertyDatabase().StorageKey("concentration")), p_key(sg.ReferencePropertyDatabase().StorageKey("porosity")); if ( first_time ) { first_time = false; // get x-coordinates at which flow is monitored std::vector<fT> xyz(4); sg.Dimensions( xyz ); fT dx = std::fabs(xyz[1]-xyz[0]); fT x2[4]; x2[0] = 0.25; x2[1] = 0.5; x2[2] = 0.75; x2[3] = 0.999; std::stringstream oss; std::string t1("BTC_x_"), t2, btc, txt(".txt"); fT a, xsec; VectorVariable<fT,dim> v; csp::Index a_key(sg.ReferencePropertyDatabase().StorageKey("aperture")), v_key(sg.ReferencePropertyDatabase().StorageKey("velocity")); // non-dimensionalised time if ( dim == 2 ) xsec = xyz[3]-xyz[2]; else xsec = (xyz[3] - xyz[2])*(xyz[5] - xyz[4]); time0 = xyz[1]-xyz[0]; time0 /= flux/(xsec*effectivePorosity(sg)); // loop over elements and identify those which have an x-coordinate above within bool largermin, smallermax; for ( stl_index i=0; i<eid.size(); i++ ) { vol[i] = 0.0; x[i] = xyz[0] + dx * x2[i]; for ( typename std::deque<Element<fT,dim> >::iterator it = sg.ReferenceMesh().ElementsBegin(); it != sg.ReferenceMesh().ElementsEnd(); it++ ) { // loop over nodes and check if x-coordinate of FE node is larger or smaller than desired x-coordinate largermin = smallermax = false; for ( stl_index j=0; j<it->Nodes(); j++ ) { if ( it->ConnectedNode(j)->x() > x[i] ) largermin = true; if ( it->ConnectedNode(j)->x() < x[i] ) smallermax = true; } // if nodes are within desired x-coordinate, backup FE and compute volume if ( largermin and smallermax ) { // compute flux across BTC plane for weighting a = it->Read( sg.ReferencePropertyStorage(), a_key.index ); it->Read( sg.ReferencePropertyStorage(), v_key.index, v ); if ( a != 0.0 ) v /= a; vel[i].push_back(v[0]); eid[i].push_back(it->ID()-1); // sum up FE area and flux vol[i] += it->Volume() * it->Read( sg.ReferencePropertyStorage(), p_key.index ) * v[0]; } } oss << x[i]; t2 = oss.str(); btc = t1+t2+txt; oss.str(""); ofs[i] = new std::ofstream( btc.c_str(), ios::out|ios::trunc ); *ofs[i] << "Time [s]\tTime [-]\tConc [-]" << endl; } } // loop over elements at x-coordinate and compute concentration in element, sum concentrations and average over total volume extern csp_float global_time; fT ctot; ScalarVariable<fT> c; typename std::list<fT>::iterator it2; for ( stl_index i=0; i<eid.size(); i++ ) { ctot = 0.0; it2 = vel[i].begin(); for ( typename std::list<stl_index>::iterator it = eid[i].begin(); it != eid[i].end(); it++, it2++ ) { // concentration weighted by area and flux across the BTC plane sg.ReferenceMesh().E(*it).PropertyValueAtBaryCenter( sg.ReferencePropertyStorage(), c_key, c ); c *= sg.ReferenceMesh().E(*it).Volume(); c *= sg.ReferenceMesh().E(*it).Read( sg.ReferencePropertyStorage(), p_key.index ); c *= *it2; ctot += c(); } *ofs[i] << global_time << "\t" << global_time/time0 << "\t" << ctot/vol[i] << endl; } } template void monitorBTC( SuperGroup<csp_float,2>& sg, csp_float flux ); template void monitorBTC( SuperGroup<csp_float,3>& sg, csp_float flux );
Example FE algorithm
Algorithm<csp_float,2> fluid_pressure; NumIntegral_dNT_op_dN_dV<csp_float,2> conductance( model.ReferencePropertyDatabase(), "conductivity", "fluid pressure", "fluid pressure" ); NumIntegral_NT_op_N_dV<csp_float,2> source( model.ReferencePropertyDatabase(), "total source", "fluid pressure" ); PointSource_rhsop<csp_float,2> boundary_flow( model.ReferencePropertyDatabase(), "influx", "fluid pressure" ); VelocityAndVolumeFlux<csp_float,2> velocity( model, "conductivity", "porosity", "fluid pressure", "fluid density", false ); fluid_pressure.Add( &conductance ); fluid_pressure.Add( &source ); fluid_pressure.Add( &boundary_flow ); fluid_pressure.AddPostProcess( &velocity );
Example function (2D) that translates a velocity applied at the LEFT boundary element by element to a flux term (velocity times area) that can be used with the PointSource PDE operator. Must be called after the initial and boundary conditions are applied and before the FE algorithm is executed.
void assignFluxToPointSource2D( SuperGroup<csp_float, 2>& sg, const char* flux ) { // backup the fluxes assigned in the configuration file in a temporary variable std::string temp_flux("temporary flux"); PropertyHandle<csp_float,2> temp_flux_var( sg, temp_flux.c_str(), SCALAR, NODE ); sg.CopyReplace( flux, temp_flux.c_str() ); // copy temporary flux (i.e., the velocity) to a temporary variable sg.InputUniformScalarValue(flux.c_str(); 0.0); // set resulting flux to zero initially csp::Index tf_key(sg.ReferencePropertyDatabase().StorageKey(temp_flux.c_str())), f_key(sg.ReferencePropertyDatabase().StorageKey(flux)); ScalarVariable<csp_float> tf, f; csp_float x[2], area, length(0.0); stl_index j; csp::String str; cout << "\nassignFluxToPointSource: Translating '" << flux << "' into a nodal point source" << endl; sg.InputUniformScalarValue( flux, 0.0 ); // set to zero initially // loop over all elements (only 2D quads and triangles) for ( std::deque<Element<csp_float, 2> >::iterator it = sg.ReferenceMesh().ElementsBegin(); it != sg.ReferenceMesh().ElementsEnd(); it++ ) { // only for elements at LEFT boundary if ( it->AtBoundary() == LEFT or it->AtBoundary() == CNR1 or it->AtBoundary() == CNR4 ) { j = 0; // first loop to calculate area for ( stl_index i=0; i<it->Nodes(); i++ ) { if ( it->ConnectedNode(i)->AtBoundary() == LEFT or it->ConnectedNode(i)->AtBoundary() == CNR1 or it->ConnectedNode(i)->AtBoundary() == CNR4 ) { x[j] = it->ConnectedNode(i)->y(); j++; } } // simple check -- quads and triangles must have two nodes at boundary if ( j != 2 ) { str = it->AtBoundary(); cerr << "\nassignFluxToPointSource: Counted less than two boundary nodes for element " << it->ID(); cerr << "\nElement at boundary: " << str.CharPointer() << endl; area = 0.0; } else { area = x[0] - x[1]; if ( area < 0.0 ) area *= -1.0; length += area; area /= static_cast<csp_float>(j); // 2 nodes per triangle or quadrilateral } // second loop to calculate nodal flux for ( stl_index i=0; i<it->Nodes(); i++ ) { if ( it->ConnectedNode(i)->AtBoundary() == LEFT or it->ConnectedNode(i)->AtBoundary() == CNR1 or it->ConnectedNode(i)->AtBoundary() == CNR4 ) { // read velocity at node i tf = it->ConnectedNode( i )->Read( sg.ReferencePropertyStorage(), tf_key.index ); // read existing result flux at node i f = it->ConnectedNode( i )->Read( sg.ReferencePropertyStorage(), f_key.index ); // add new flux to this value f += tf() * area; // store it back to node it->ConnectedNode( i )->Store( sg.ReferencePropertyStorage(), f_key.index, f ); } } } } printRangeOfVariable( sg, flux ); cout << "\nassignFluxToPointSource: Successfully assigned '" << flux << "' to the nodes... " << endl; }
Here are some general C++ code snippets:
Well-known C functions like "atoi" are deprecated in C++ for the following reason. When "atoi" encounters a string with no numerical sequence, it returns zero (0). If the string holds a valid sequence of digits that represents the number 0, it also returns a 0, making it impossible to tell from the return value alone whether the string holds a valid number or not. A valid C++ code that utilizes STL and templates and avoids the described drawback is presented below:
template <class T> bool from_string(T& t, const std::string& s, std::ios_base& (*f)(std::ios_base&)) { std::istringstream iss(s); return !(iss >> f >> t).fail(); }
This single function allows one to parse a string into any type of number representation due to the power of templates. The last argument makes it possible to treat a string not only as a decimal number, but also as hexadecimal, octal, etc.
The example of its usage looks like this:
int D; string str = "FF"; from_string<int>(D, str, std::hex); cout << D; ------ 255
You can, of course, use the usual for-loop to do this, but somewhat fancier one-liner you can try instead is this:
#include <iterator> .... copy (vect.begin(), vect.end(), ostream_iterator<int> (cout, " "));
Here "vect" is a vector<int> and the very last string argument is a separator that will be printed in between the vector elements.
