Select Git revision
LineEdgeBasis.cpp
Forked from
gmsh / gmsh
Source project has a limited visibility.
-
Nicolas Marsic authored
Reorganisation of FunctionSpace -- Basis can be ONLY evaluated + Derivatives -- New Dof ordering: Iterate on Entity and then on Function Order
Nicolas Marsic authoredReorganisation of FunctionSpace -- Basis can be ONLY evaluated + Derivatives -- New Dof ordering: Iterate on Entity and then on Function Order
HighOrder.cpp 60.03 KiB
// Gmsh - Copyright (C) 1997-2019 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
//
// Contributor(s):
// Koen Hillewaert
//
#include <sstream>
#include <vector>
#include "GmshConfig.h"
#include "HighOrder.h"
#include "MLine.h"
#include "MTriangle.h"
#include "MQuadrangle.h"
#include "MTetrahedron.h"
#include "MHexahedron.h"
#include "MPrism.h"
#include "MPyramid.h"
#include "GmshMessage.h"
#include "OS.h"
#include "fullMatrix.h"
#include "BasisFactory.h"
#include "InnerVertexPlacement.h"
#include "Context.h"
#if defined(HAVE_OPTHOM)
#include "HighOrderMeshFastCurving.h"
#include "HighOrderMeshPeriodicity.h"
#endif
// Functions that help optimizing placement of points on geometry
// The aim here is to build a polynomial representation that consist
// in polynomial segments of equal length
static double mylength(GEdge *ge, int i, double *u)
{
return ge->length(u[i], u[i + 1], 10);
}
static void myresid(int N, GEdge *ge, double *u, fullVector<double> &r,
double *weight = NULL)
{
double L[100];
for(int i = 0; i < N - 1; i++) L[i] = mylength(ge, i, u);
if(weight)
for(int i = 0; i < N - 2; i++)
r(i) = L[i + 1] / weight[i + 1] - L[i] / weight[i];
else
for(int i = 0; i < N - 2; i++) r(i) = L[i + 1] - L[i];
}
static bool computeEquidistantParameters(GEdge *ge, double u0, double uN, int N,
double *u, double underRelax)
{
const double PRECISION = 1.e-6;
const int MAX_ITER = 50;
const double eps = (uN - u0) * 1.e-5;
// newton algorithm
// N is the total number of points (3 for quadratic, 4 for cubic ...)
// u[0] = u0;
// u[N-1] = uN;
// initialize as equidistant in parameter space
u[0] = u0;
double du = (uN - u0) / (N - 1);
for(int i = 1; i < N; i++) {
u[i] = u[i - 1] + du;