// Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced // at the Lawrence Livermore National Laboratory. All Rights reserved. See files // LICENSE and NOTICE for details. LLNL-CODE-806117. // // This file is part of the MFEM library. For more information and source code // availability visit https://mfem.org. // // MFEM is free software; you can redistribute it and/or modify it under the // terms of the BSD-3 license. We welcome feedback and contributions, see file // CONTRIBUTING.md for details. // // ------------------------------------------------------------ // NURBS CurveInt Miniapp: Interpolate a Curve in a NURBS Patch // ------------------------------------------------------------ // // Compile with: make nurbs_curveint // // Sample runs: ./nurbs_curveint -uw -n 9 // ./nurbs_curveint -nw -n 9 // // Description: This example code demonstrates the use of MFEM to interpolate a // curve in a NURBS patch. We first define a square shaped NURBS // patch. We then interpolate a sine function on the bottom // edge. The results can be viewed in VisIt. // // We use curve interpolation for curves with all weights being 1, // B-splines, and curves with not all weights being 1, NURBS. The // spacing in both cases is chosen differently. #include "mfem.hpp" #include #include using namespace std; using namespace mfem; KnotVector *UniformKnotVector(int order, int ncp) { if (order>=ncp) { mfem_error("UniformKnotVector: ncp should be at least order + 1"); } KnotVector *kv = new KnotVector(order, ncp); for (int i = 0; i < order+1; i++) { (*kv)[i] = 0.0; } for (int i = order+1; i < ncp; i++) { (*kv)[i] = (i-order)/double(ncp-order); } for (int i = ncp ; i < ncp + order + 1; i++) { (*kv)[i] = 1.0; } return kv; } int main(int argc, char *argv[]) { // Parse command-line options. OptionsParser args(argc, argv); double l = 1.0; double a = 0.1; int ncp = 9; int order = 2; bool ifbspline = true; bool visualization = true; bool visit = true; args.AddOption(&l, "-l", "--box-side-length", "Height and width of the box"); args.AddOption(&a, "-a", "--sine-ampl", "Amplitude of the fitted sine function."); args.AddOption(&ncp, "-n", "--ncp", "Number of control points used over four box sides."); args.AddOption(&order, "-o", "--order", "Order of the NURBSPatch"); args.AddOption(&ifbspline, "-uw", "--unit-weight", "-nw", "--non-unit-weight", "Use a unit-weight for B-splines (default) or not: for general NURBS"); args.AddOption(&visualization, "-vis", "--visualization", "-no-vis", "--no-visualization", "Enable or disable GLVis visualization. This is a dummy option to enable testing."); args.AddOption(&visit, "-visit", "--visit", "-no-visit", "--no-visit", "Enable or disable VisIt visualization."); // Parse and print command line options args.Parse(); if (!args.Good()) { args.PrintUsage(cout); return 1; } args.PrintOptions(cout); if (order < 2 && ifbspline == false) { mfem_error("For a non unity weight, the order should be at least 2."); } KnotVector *kv_o1 = UniformKnotVector(1, 2); KnotVector *kv = UniformKnotVector(order, ncp); // 1. Create a box shaped NURBS patch NURBSPatch patch(kv_o1, kv_o1, 3); // Set weights for (int j = 0; j < 2; j++) for (int i = 0; i < 2; i++) { patch(i,j,2) = 1.0; } // Define patch corners which are box corners patch(0,0,0) = -0.5*l; patch(0,0,1) = -0.5*l; patch(1,0,0) = 0.5*l; patch(1,0,1) = -0.5*l; patch(0,1,0) = -0.5*l; patch(0,1,1) = 0.5*l; patch(1,1,0) = 0.5*l; patch(1,1,1) = 0.5*l; // 2. Interpolation process Array xy(2); xy[0] = new Vector(); xy[1] = new Vector(); Vector xi_args, u_args; Array i_args; xy[0]->SetSize(ncp); xy[1]->SetSize(ncp); // Refine direction which has fitting if (!ifbspline) { // We alter the weight for demonstration purposes to a random value. This // is not necessary for general curve fitting. patch.DegreeElevate(0, 1); patch(1,0,2) = sqrt(2)/2; patch.DegreeElevate(0, order-kv_o1->GetOrder()-1); } else { patch.DegreeElevate(0, order-kv_o1->GetOrder()); } patch.KnotInsert(0, *kv); // We locate the control points at the location of the maxima of the // knot vectors. This works very well for patches with unit weights. kv->FindMaxima(i_args,xi_args, u_args); for (int i = 0; i < ncp; i++) { (*xy[0])[i] = u_args[i]*l; (*xy[1])[i] = a * sin((*xy[0])[i]/l*2*M_PI)-0.5*l; (*xy[0])[i] -= 0.5*l; } kv->FindInterpolant(xy); // Apply interpolation to patch for (int i = 0; i < ncp; i++) { patch(i,0,0) = (*xy[0])[i]; patch(i,0,1) = (*xy[1])[i]; } if (!ifbspline) { // Convert to homogeneous coordinates. FindInterpolant returns // Cartesian coordinates. for (int i = 0; i < ncp; i++) { patch(i,0,0) *= patch(i,0,2); patch(i,0,1) *= patch(i,0,2); } } // Refinement in curve interpolation direction patch.DegreeElevate(1, order-kv_o1->GetOrder()); patch.KnotInsert(1, *kv); // 3. Open and write mesh output file string mesh_file("sin-fit.mesh"); ofstream output(mesh_file.c_str()); output<<"MFEM NURBS mesh v1.0"<PrintInfo(); if (visit) { // Print mesh to file for visualization VisItDataCollection dc = VisItDataCollection("mesh", mesh); dc.SetPrefixPath("CurveInt"); dc.SetCycle(0); dc.SetTime(0.0); dc.Save(); } delete mesh; delete kv_o1; delete kv; delete xy[0]; delete xy[1]; return 0; }