SmoothSurfSet.h

//===========================================================================
// GoTools - SINTEF Geometry Tools version 1.1
//
// GoTools module: CORE
//
// Copyright (C) 2000-2007 SINTEF ICT, Applied Mathematics, Norway.
//
// This program is free software; you can redistribute it and/or          
// modify it under the terms of the GNU General Public License            
// as published by the Free Software Foundation version 2 of the License. 
//
// This program is distributed in the hope that it will be useful,        
// but WITHOUT ANY WARRANTY; without even the implied warranty of         
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          
// GNU General Public License for more details.                           
//
// You should have received a copy of the GNU General Public License      
// along with this program; if not, write to the Free Software            
// Foundation, Inc.,                                                      
// 59 Temple Place - Suite 330,                                           
// Boston, MA  02111-1307, USA.                                           
//
// Contact information: E-mail: tor.dokken@sintef.no                      
// SINTEF ICT, Department of Applied Mathematics,                         
// P.O. Box 124 Blindern,                                                 
// 0314 Oslo, Norway.                                                     
//
// Other licenses are also available for this software, notably licenses
// for:
// - Building commercial software.                                        
// - Building software whose source code you wish to keep private.        
//===========================================================================
#ifndef _GOSMOOTHSURFSET_H_
#define _GOSMOOTHSURFSET_H_


//   -----------------------------------------------------------------------
//      Interface file for class SmoothSurfSET
//   -----------------------------------------------------------------------
//
//       Used to modify a tensor product B-spline surface with respect
//       to conditions on smoothness, editing constraints and boundary
//       conditions.
//
//       Implementation of the member functions are given in the
//       following files:
//
//          1. SmoothSurfSet.C
//
//   -----------------------------------------------------------------------
//    Written by: Vibeke Skytt                            04.2002.
//   -----------------------------------------------------------------------

#include "SplineSurface.h"
#include "ConstraintDefinitions.h"
#include <vector>

namespace Go
{


class SmoothSurfSet
{


public:
    SmoothSurfSet();

    SmoothSurfSet(bool copy_coefs);

    virtual
    ~SmoothSurfSet();

    void attach(std::vector<boost::shared_ptr<SplineSurface> >& insf, 
                std::vector<std::vector<int> >& coef_known,
                int num_side_constraints = 0, 
                int has_normal_cond = 0);

    // @@@ VSK. Can it be relevant to use different weights for 
    // different surfaces, and how to define the weights in that case?
    virtual
    void setOptimize(double weight1, double weight2, double weight3);

    void setLeastSquares(std::vector<std::vector<double> >& pnts,
                         std::vector<std::vector<double> >& param_pnts,
                         std::vector<std::vector<double> >& pnt_weights,
                         double weight);

    int setNormalCond(std::vector<std::vector<double> >& pnts,
                      std::vector<std::vector<double> >& param_pnts,
                      std::vector<std::vector<double> >& pnt_weights,
                      double weight);

    void approxOrig(double weight);

    void setSideConstraints(std::vector<sideConstraintSet>& constraints);

    // We may have side constraints which are not suitable for exact equality as
    // spline solution space may not be large enough. We therefore allow using least
    // squares to minimize the error.
    // This applies in particular to constraint involving higher order derivatives.
    // @@ Currently function uses newmat. As matrix is sparce we may do it another way.
    void setApproxSideConstraints(std::vector<sideConstraintSet>& constraints,
                                  double weight);

    int equationSolve(std::vector<boost::shared_ptr<SplineSurface> >& surfaces);


protected:
    int kdim_;             // 3 if normal conditions are included, otherwise 1.
    int idim_;             // Dimension of geometry space. 
    int idim1_;            // Dimension of projective space.
    int ider_;             // Maximum derivative involved in the computations. 

    int kncond_;           // Number of unknown coefficents + number of constraints.
    int knconstraint_;     // Number of side constraints.

    std::vector<std::vector<int>::iterator> coef_known_;
    // Array indicating the status of coefficients, i.e.
    // free, fixed, not involved, equal to a given
    // coefficients.
    std::vector<std::vector<int> > pivot_; // Array giving the position of the
    // free coefficients in the equation system.

    // Parameters defining the spline space.
    std::vector<boost::shared_ptr<SplineSurface> > srfs_;   // Pointer to input surfaces.

    const int copy_coefs_; // Whether we are to copy the coefs from input surfaces
                           // (or work directly on coefficients in input surfaces).

    // Parameters used to define the specific input spline surfaces.
    std::vector<std::vector<double> > coef_array_;

    //    std::vector<double>::iterator  scoef;   // Pointer to surface coefficients.      

    // Storage of the equation system.
    std::vector<double> gmat_;       // Matrix at left side of equation system.
    std::vector<double> gright_;     // Right side of equation system.      

    virtual
    void getBasis(double *sb1, double *sb2, int kk1, int kk2, 
                  int kleft1, int kleft2, int ider, double *sbasis);


private:

    typedef struct integralInfo
    {
        // Parameters used in integration
        std::vector<double> vec1, vec2;
        double ***integral1;  // Array used to store integrals of inner product
        // of derivatives of B-splines in 1. par. dir.
        // The 1st index runs over the derivatives, the 2nd & 3rd in u- and v-dir.
        double ***integral2;  // Array used to store integrals of inner product
        // of derivatives of B-splines in 2. par. dir.
        // The 1st index runs over the derivatives, the 2nd & 3rd in u- and v-dir.
        bool integralset; // Whether integral1 & integral2 have been computed.
        int der; // The number of derivatives to compute.
      
        integralInfo() 
        { integral1 = 0; integral2 = 0; integralset = false; der = -1; }

        ~integralInfo()
        { erase(); }

        void resize(int ider, int in1, int in2)
        {
            int ki, kj;
            vec1.resize((ider+1)*in1*in1);
            vec2.resize((ider+1)*in2*in2);
            std::fill(vec1.begin(), vec1.end(), 0.0);
            std::fill(vec2.begin(), vec2.end(), 0.0);

            integral1 = new double**[ider+1];
            integral2 = new double**[ider+1];

            for (ki=0; ki<=ider; ki++)
                {
                    integral1[ki] = new double*[in1];
                    integral2[ki] = new double*[in2];

                    for (kj=0; kj<in1; kj++)
                        integral1[ki][kj] = &vec1[(ki*in1+kj)*in1];

                    for (kj=0; kj<in2; kj++)
                        integral2[ki][kj] = &vec2[(ki*in2+kj)*in2];
                }

            integralset = 0;
            der = ider;
        }

        void erase()
        {
            int ki;
            for (ki=0; ki<=der; ki++)
                {
                    delete [] integral1[ki];
                    delete [] integral2[ki];
                }
            delete [] integral1;
            delete [] integral2;
            integralset = false;
        }
    } integralInfo;

    std::vector<integralInfo> surf_integral_; // Integral calculations for the surfaces.

    void removeKnownCoefs(std::vector<sideConstraintSet>& constraints) const;

};

}

#endif //

---