CarveMeshPtr create_mesh(const MatrixFr& vertices, const MatrixIr& faces) {
        const size_t num_vertices = vertices.rows();
        const size_t num_faces = faces.rows();

        if (vertices.cols() != 3) {
            throw NotImplementedError("Only 3D mesh is supported.");
        }
        if (faces.cols() != 3) {
            throw NotImplementedError("Only triangle mesh is supported.");
        }

        std::vector<CarveVector> points;

        for (size_t i=0; i<num_vertices; i++) {
            const auto& v = vertices.row(i);
            CarveVector p;
            p.v[0] = v[0];
            p.v[1] = v[1];
            p.v[2] = v[2];
            points.push_back(p);
        }

        std::vector<int> raw_faces;
        raw_faces.reserve(num_faces * 4);
        for (size_t i=0; i<num_faces; i++) {
            raw_faces.push_back(3);
            raw_faces.push_back(faces(i,0));
            raw_faces.push_back(faces(i,1));
            raw_faces.push_back(faces(i,2));
        }

        return CarveMeshPtr(new CarveMesh(points, num_faces, raw_faces));
    }

FESettingFactory& FESettingFactory::with_material(
        const std::string& material_name) {
    size_t dim = m_mesh->getDim();
    if (material_name == "test_material") {
        m_material = Material::create_isotropic(dim, 1.0, 1.0, 0.0);
    } else if (material_name == "periodic_material") {
        MaterialPtr mat1 = Material::create_isotropic(dim, 1.0, 1.0, 0.0);
        MaterialPtr mat2 = Material::create_isotropic(dim, 1.0, 2.0, 0.0);
        VectorF axis = VectorF::Zero(dim);
        axis[0] = 1;
        m_material = Material::create_periodic(mat1, mat2, axis, 1.0, 0.5, 0.0);
    } else if (material_name == "homogenized_material") {
        if (dim == 3) {
            throw NotImplementedError("Homogenized material is not supported in 3D");
        }
        size_t tensor_size = dim * dim;
        MatrixF C(tensor_size, tensor_size);
        C << 4.0/3.0,     0.0,     0.0, 7.0/5.0,
             0.0    ,     0.0,     0.0,     0.0,
             0.0    ,     0.0,     0.0,     0.0,
             7.0/5.0,     0.0,     0.0, 3.0/2.0;
        m_material = Material::create(1.0, C);
    } else {
        std::stringstream err_msg;
        err_msg << "Material \"" << material_name << "\" is not supported.";
        throw NotImplementedError(err_msg.str());
    }
    return *this;
}

Elements::Ptr Elements::adapt_boundary(Mesh::Ptr mesh) {
    if (mesh->get_num_voxels() > 0) {
        const size_t vertex_per_voxel = mesh->get_vertex_per_voxel();
        switch (vertex_per_voxel) {
            case 4:
                return Ptr(new TriangleElements(mesh));
            default:
                std::stringstream err_msg;
                err_msg << "Voxel with " << vertex_per_voxel
                    << " vertices is not supported yet.";
                throw NotImplementedError(err_msg.str());
        }
    } else {
        const size_t vertex_per_face = mesh->get_vertex_per_face();
        switch (vertex_per_face) {
            case 3:
                return Ptr(new EdgeElements(mesh));
            default:
                std::stringstream err_msg;
                err_msg << "Face with " << vertex_per_face
                    << " vertices is not supported yet.";
                throw NotImplementedError(err_msg.str());
        }
    }
}

void PointLocator::init_elements() {
    const size_t dim = m_mesh->get_dim();

    if (dim == 2) {
        if (m_mesh->get_num_faces() == 0) {
            throw RuntimeError("2D Mesh has no faces.");
        }
        m_elements = m_mesh->get_faces();
        m_vertex_per_element = m_mesh->get_vertex_per_face();
        if (m_vertex_per_element != 3) {
            throw NotImplementedError(
                    "Only triangle elements are supported in 2D");
        }
    } else if (dim == 3) {
        if (m_mesh->get_num_voxels() == 0) {
            throw RuntimeError("3D Mesh has no voxels.");
        }
        m_elements = m_mesh->get_voxels();
        m_vertex_per_element = m_mesh->get_vertex_per_voxel();
        if (m_vertex_per_element != 4) {
            throw NotImplementedError(
                    "Only tetrahedron elements are supported in 2D");
        }
    } else {
        throw NotImplementedError("Only 2D and 3D mesh are supported");
    }
}

void HelpWriterContext::writeOptionItem(const std::string &name,
                                        const std::string &args,
                                        const std::string &description) const
{
    File &file = outputFile();
    switch (outputFormat())
    {
        case eHelpOutputFormat_Console:
            // TODO: Generalize this when there is need for it; the current,
            // special implementation is in CommandLineHelpWriter.
            GMX_THROW(NotImplementedError("Option item formatting for console output not implemented"));
            break;
        case eHelpOutputFormat_Man:
            file.writeLine(formatString(".BI \"\\%s\" \" %s\"", name.c_str(), args.c_str()));
            file.writeString("    ");
            writeTextBlock(description);
            file.writeLine();
            break;
        case eHelpOutputFormat_Html:
        {
            std::string substArgs =
                substituteMarkupAndWrapToString(TextLineWrapperSettings(), args);
            file.writeLine(formatString("<dt><b><tt>%s</tt></b> %s</dt>", name.c_str(),
                                        substArgs.c_str()));
            file.writeLine("<dd>");
            writeTextBlock(description);
            file.writeLine("</dd>");
            break;
        }
        default:
            GMX_THROW(NotImplementedError(
                              "This output format is not implemented"));
    }
}

RCP<const Basic> log(const RCP<const Basic> &arg)
{
    if (eq(*arg, *zero)) {
        throw NotImplementedError(
            "log(0) is complex infinity. Yet to be implemented");
    }
    if (eq(*arg, *one))
        return zero;
    if (eq(*arg, *E))
        return one;
    if (is_a_Number(*arg)) {
        RCP<const Number> _arg = rcp_static_cast<const Number>(arg);
        if (not _arg->is_exact()) {
            return _arg->get_eval().log(*_arg);
        } else if (_arg->is_negative()) {
            throw NotImplementedError(
                "Imaginary Result. Yet to be implemented");
        }
    }
    if (is_a<Rational>(*arg)) {
        RCP<const Integer> num, den;
        get_num_den(static_cast<const Rational &>(*arg), outArg(num),
                    outArg(den));
        return sub(log(num), log(den));
    }
    return make_rcp<const Log>(arg);
}

void ExcessTerm::construct(const std::vector<CoolPropFluid*> &components)
{
    std::string _model;
    N = components.size();

    F.resize(N, std::vector<double>(N, 0));
    DepartureFunctionMatrix.resize(N);

    for (unsigned int i = 0; i < N; ++i)
    {
        DepartureFunctionMatrix[i].resize(N);
        for (unsigned int j = 0; j < N; ++j)
        {
            if (i == j){ continue; }

            std::string CAS1 = components[i]->CAS;
            std::vector<std::string> CAS(2,"");
            CAS[0] = components[i]->CAS;
            CAS[1] = components[j]->CAS;
            std::sort(CAS.begin(), CAS.end());

            std::vector<Dictionary> & vd = mixtureexcesslibrary.excess_map[CAS];
            if (vd.size() != 1) { throw NotImplementedError(); }
            // Get a reference to the dictionary itself to save a few dereferences
            Dictionary &dic = vd[0];

            std::string model = dic.get_string("model");

            if (!model.compare("Kunz-JCED-2012"))
            {
                F[i][j] = dic.get_number("F");

                std::vector<double> n = dic.get_double_vector("n");
                std::vector<double> d = dic.get_double_vector("d");
                std::vector<double> t = dic.get_double_vector("t");
                // Terms for the gaussian
                std::vector<double> eta = dic.get_double_vector("eta");
                std::vector<double> epsilon = dic.get_double_vector("epsilon");
                std::vector<double> beta = dic.get_double_vector("beta");
                std::vector<double> gamma = dic.get_double_vector("gamma");
                int Npower = static_cast<int>(dic.get_number("Npower"));
                DepartureFunctionMatrix[i][j].reset(new GERG2008DepartureFunction(n,d,t,eta,epsilon,beta,gamma,Npower));
            }
            else if (!model.compare("Lemmon-JPCRD-2004") || !model.compare("Lemmon-JPCRD-2000"))
            {
                throw NotImplementedError();
            }
            else
            {
                throw ValueError();
            }
        }
    }
}

    std::vector<VectorI> enumerate(const VectorI& repetitions) {
        std::vector<VectorI> result;
        const size_t dim = repetitions.size();
        if (dim == 2) {
            for (size_t i=0; i<repetitions[0]; i++) {
                for (size_t j=0; j<repetitions[1]; j++) {
                    result.push_back(Vector2I(i,j));
                }
            }
        } else if (dim == 3) {
            for (size_t i=0; i<repetitions[0]; i++) {
                for (size_t j=0; j<repetitions[1]; j++) {
                    for (size_t k=0; k<repetitions[2]; k++) {
                        result.push_back(Vector3I(i,j,k));
                    }
                }
            }
        } else {
            std::stringstream err_msg;
            err_msg << "Unsupported dim: " << dim;
            throw NotImplementedError(err_msg.str());
        }

        return result;
    }

/// Use the single-phase table to evaluate an output
double CoolProp::TTSEBackend::evaluate_single_phase_derivative(SinglePhaseGriddedTableData &table, parameters output, double x, double y, std::size_t i, std::size_t j, std::size_t Nx, std::size_t Ny)
{
	if (Nx == 1 && Ny == 0){
		if (output == table.xkey) { return 1.0; }
		if (output == table.ykey) { return 0.0; }
	}
	else if (Ny == 1 && Nx == 0){
		if (output == table.ykey) { return 1.0; }
		if (output == table.xkey) { return 0.0; }
	}
    
    connect_pointers(output, table);
    
    // Distances from the node
	double deltax = x - table.xvec[i];
    double deltay = y - table.yvec[j];
    double val;
    // Calculate the output value desired
    if (Nx == 1 && Ny == 0){
		if (output == table.xkey) { return 1.0; }
		if (output == table.ykey) { return 0.0; }
		val = (*dzdx)[i][j] + deltax*(*d2zdx2)[i][j] + deltay*(*d2zdxdy)[i][j];
	}
	else if (Ny == 1 && Nx == 0){
		if (output == table.ykey) { return 1.0; }
		if (output == table.xkey) { return 0.0; }
		val = (*dzdy)[i][j] + deltay*(*d2zdy2)[i][j] + deltax*(*d2zdxdy)[i][j];
	}
	else{
		throw NotImplementedError("only first derivatives currently supported");
	}
    return val;
}

SparseSolverFactory::SparseSolverFactory(const std::string& solver_type) {
    if (solver_type == "LDLT") {
        m_solver = SparseSolverPtr(
                new SparseSolverImplementation<Eigen::SimplicialLDLT<ZSparseMatrix::ParentType> >);
    } else if (solver_type == "LLT") {
        m_solver = SparseSolverPtr(
                new SparseSolverImplementation<Eigen::SimplicialLLT<ZSparseMatrix::ParentType> >);
    } else if (solver_type == "CG") {
        m_solver = SparseSolverPtr(
                new SparseSolverImplementation<Eigen::ConjugateGradient<ZSparseMatrix::ParentType> >);
    } else if (solver_type == "SparseLU") {
        m_solver = SparseSolverPtr(
                new SparseSolverImplementation<Eigen::SparseLU<ZSparseMatrix::ParentType> >);
    } else if (solver_type == "UmfPackLU") {
        m_solver = SparseSolverPtr(
                new SparseSolverImplementation<Eigen::UmfPackLU<ZSparseMatrix::ParentType> >);
    } else if (solver_type == "UmfPack") {
        m_solver = SparseSolverPtr(
                new SparseSolverImplementation<UmfpackFactorizer>);
    } else {
        std::stringstream err_msg;
        err_msg << "Unsupported solver type " << solver_type;
        throw NotImplementedError(err_msg.str());
    }
}

std::unique_ptr<BoxDeformation>
prepareBoxDeformation(const matrix     &initialBox,
                      t_commrec        *cr,
                      const t_inputrec &inputrec)
{
    if (!inputrecDeform(&inputrec))
    {
        return nullptr;
    }
    if (!EI_DYNAMICS(inputrec.eI))
    {
        GMX_THROW(NotImplementedError("Box deformation is only supported with dynamical integrators"));
    }

    matrix box;
    // Only the rank that read the tpr has the global state, and thus
    // the initial box, so we pass that around.
    if (SIMMASTER(cr))
    {
        copy_mat(initialBox, box);
    }
    if (PAR(cr))
    {
        gmx_bcast(sizeof(box), box, cr);
    }

    return compat::make_unique<BoxDeformation>(inputrec.delta_t,
                                               inputrec.init_step,
                                               inputrec.deform,
                                               box);
}

/** If the fluid type is mole-based, it does not do anything. Otherwise,
 *  it converts the mole fraction to the required input. */
double IncompressibleFluid::inputFromMole (double T,     double x){
    if (this->xid==IFRAC_PURE) {
            return _HUGE;
    } else if (this->xid==IFRAC_MOLE) {
        return x;
    } else {
        throw NotImplementedError("Mole composition conversion has not been implemented.");
        /*
        switch (mole2input.type) {
            case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
                return poly.evaluate(mole2input.coeffs, T, x, 0, 0, 0.0, 0.0); // TODO: make sure Tbase and xbase are defined in the correct way
                break;
            case IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL:
                return baseExponential(mole2input, x, 0.0);
                break;
            case IncompressibleData::INCOMPRESSIBLE_LOGEXPONENTIAL:
                return baseLogexponential(mole2input, x, 0.0);
                break;
            case IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL:
                return exp(poly.evaluate(mole2input.coeffs, T, x, 0, 0, 0.0, 0.0)); // TODO: make sure Tbase and xbase are defined in the correct way
                break;
            case IncompressibleData::INCOMPRESSIBLE_POLYOFFSET:
                return basePolyOffset(mole2input, T, x);
                break;
            case IncompressibleData::INCOMPRESSIBLE_NOT_SET:
                throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,mole2input.type));
                break;
            default:
                throw ValueError(format("%s (%d): Your function type \"[%d]\" is unknown.",__FILE__,__LINE__,mole2input.type));
                break;
        }
        return _HUGE;
        */
    }
}

VectorF extract_vertices(std::unique_ptr<DracoMesh>& draco_mesh, size_t& dim) {
    const auto num_vertices = draco_mesh->num_points();
    const auto positions = draco_mesh->GetNamedAttribute(
            draco::GeometryAttribute::POSITION);
    assert(positions->IsValid());
    dim = positions->num_components();

    VectorF vertices(num_vertices * dim);
    if (dim == 2) {
        for (size_t i=0; i<num_vertices; i++) {
            const auto p = positions->template GetValue<Float, 2>(
                    draco::AttributeValueIndex(i));
            vertices(i*2  ) = p[0];
            vertices(i*2+1) = p[1];
        }
    } else if (dim == 3) {
        for (size_t i=0; i<num_vertices; i++) {
            const auto p = positions->template GetValue<Float, 3>(
                    draco::AttributeValueIndex(i));
            vertices(i*3  ) = p[0];
            vertices(i*3+1) = p[1];
            vertices(i*3+2) = p[2];
        }
    } else {
        throw NotImplementedError("Draco mesh encodes high dimensional data");
    }
    return vertices;
}

void HelpWriterContext::writeOptionItem(const std::string &name,
                                        const std::string &value,
                                        const std::string &defaultValue,
                                        const std::string &info,
                                        const std::string &description) const
{
    TextWriter &file = outputFile();
    switch (outputFormat())
    {
        case eHelpOutputFormat_Console:
        {
            TextTableFormatter &formatter(impl_->state_->consoleOptionsFormatter());
            formatter.clear();
            formatter.addColumnLine(0, name);
            formatter.addColumnLine(1, value);
            if (!defaultValue.empty())
            {
                formatter.addColumnLine(2, "(" + defaultValue + ")");
            }
            if (!info.empty())
            {
                formatter.addColumnLine(3, "(" + info + ")");
            }
            TextLineWrapperSettings settings;
            settings.setIndent(11);
            settings.setLineLength(78);
            std::string formattedDescription
                = substituteMarkupAndWrapToString(settings, description);
            file.writeLine(formatter.formatRow());
            file.writeLine(formattedDescription);
            break;
        }
        case eHelpOutputFormat_Rst:
        {
            std::string args(value);
            if (!defaultValue.empty())
            {
                args.append(" (");
                args.append(defaultValue);
                args.append(")");
            }
            if (!info.empty())
            {
                args.append(" (");
                args.append(info);
                args.append(")");
            }
            file.writeLine(formatString("``%s`` %s", name.c_str(), args.c_str()));
            TextLineWrapperSettings settings;
            settings.setIndent(4);
            file.writeLine(substituteMarkupAndWrapToString(settings, description));
            break;
        }
        default:
            GMX_THROW(NotImplementedError(
                              "This output format is not implemented"));
    }
}

OffsetParameters::OffsetParameters(WireNetwork::Ptr wire_network,
        OffsetParameters::TargetType type,
        Float default_value) :
    m_wire_network(wire_network),
    m_type(type),
    m_default_offset(default_value) {
        if (m_type == ParameterCommon::EDGE) {
            throw NotImplementedError("Only vertex offset is supported.");
        }
    }

void MeshGeometry::extract_faces_from_voxels() {
    if (m_vertex_per_voxel == 4) {
        extract_faces_from_tets();
    } else if (m_vertex_per_voxel == 8) {
        extract_faces_from_hexes();
    } else {
        std::stringstream err_msg;
        err_msg << "Unsupported voxel type with " << m_vertex_per_voxel << " vertices per voxel";
        throw NotImplementedError(err_msg.str());
    }
}

HelpWriterContext::HelpWriterContext(File *file, HelpOutputFormat format)
    : impl_(new Impl(file, format))
{
    if (format != eHelpOutputFormat_Console)
    {
        // TODO: Implement once the situation with Redmine issue #969 is more
        // clear.
        GMX_THROW(NotImplementedError(
                    "This output format is not implemented"));
    }
}

void EigenSolver::compute_batch_symmetric(size_t dim, VectorF matrices) {
    if (dim == 2) {
        compute_batch_symmetric_2x2(matrices);
    } else if (dim == 3) {
        compute_batch_symmetric_3x3(matrices);
    } else {
        std::stringstream err_msg;
        err_msg << "EigenSolver does not support " << dim << "D matrix yet.";
        throw NotImplementedError(err_msg.str());
    }
}

void IsotropicTransforms::initialize() {
    if (m_dim == 2) {
        initialize_2D_reflections();
        initialize_2D_rotations();
    } else if (m_dim == 3) {
        initialize_3D_reflections();
        initialize_3D_rotations();
    } else {
        throw NotImplementedError("Only 2D and 3D is supported");
    }
}

UniformMaterial::UniformMaterial(
        Float density, const MatrixF& material_tensor) :
        m_density(density), m_material_tensor(material_tensor) {
    const size_t rows = m_material_tensor.rows();
    if (rows == 9) {
        m_dim = 3;
    } else if (rows == 4) {
        m_dim = 2;
    } else {
        throw NotImplementedError("Material tensor size not supported.");
    }
}