Matlab中的非线性固体力学求解器_MATLAB_下载.zip

classdef NurbsTools < handle properties(Access = private) basis_function_ nurbs_data_ end methods % Constructor function this = NurbsTools(basis_function) this.basis_function_ = basis_function; this.nurbs_data_ = basis_function.topology_data_.domain_patch_data_.nurbs_data_; end % Plot parametric mesh formed by uniqued knot vectors function plotParametricMesh(this) switch this.nurbs_data_.getGeometryDimension() case 1 unique_knot = unique(this.nurbs_data_.knot_vectors_{1}); plot(unique_knot, unique_knot*0, 'ko'); xlabel('\xi'); case 2 unique_knot_1 = unique(this.nurbs_data_.knot_vectors_{1}); unique_knot_2 = unique(this.nurbs_data_.knot_vectors_{2}); unique_knot_3 = 0; this.plotMeshPatch(unique_knot_1, unique_knot_2, unique_knot_3); xlabel('\xi'); ylabel('\eta'); case 3 unique_knot_1 = unique(this.nurbs_data_.knot_vectors_{1}); unique_knot_2 = unique(this.nurbs_data_.knot_vectors_{2}); unique_knot_3 = unique(this.nurbs_data_.knot_vectors_{3}); % Plot face zeta = 0 this.plotMeshPatch(unique_knot_1, unique_knot_2, 0); % Plot face zeta = 1 this.plotMeshPatch(unique_knot_1, unique_knot_2, 1); % Plot face eta = 0 this.plotMeshPatch(unique_knot_1, 0, unique_knot_3); % Plot face eta = 1 this.plotMeshPatch(unique_knot_1, 1, unique_knot_3); % Plot face xi = 0 this.plotMeshPatch(0, unique_knot_2, unique_knot_3); % Plot face xi = 1 this.plotMeshPatch(1, unique_knot_2, unique_knot_3); xlabel('\xi'); ylabel('\eta'); zlabel('\zeta'); end end % Evaluate nurbs or its derivatives or its hassian matrix at xi function [position, gradient, hassian] = evaluateNurbs(this, xi, varargin) if ~isequal(size(xi,2), this.nurbs_data_.getGeometryDimension()) disp('Input parametric coordinates dimension are mismatched!'); return; end import BasisFunction.IGA.QueryUnit query_unit = QueryUnit(); domain_patch = this.basis_function_.topology_data_.domain_patch_data_; % Current shape funciton have to query point by point geo_dim = this.nurbs_data_.getGeometryDimension(); position = zeros(size(xi,1), 3); gradient = cell(1,geo_dim); hassian = cell(1,geo_dim*(geo_dim+1)/2); if isempty(varargin) % default option varargin = {'position', 'gradient', 'hassian'}; end if any(strcmp(varargin, 'hassian')) for i = 1:size(xi,1) query_unit.query_protocol_ = {domain_patch, xi(i,:), 2}; this.basis_function_.query(query_unit); non_zero_id = query_unit.non_zero_id_; R = query_unit.evaluate_basis_{1}; position(i,:) = R * this.nurbs_data_.control_points_(non_zero_id,1:3); for dim_i = 1:geo_dim temp = query_unit.evaluate_basis_{2}(dim_i,:); gradient{dim_i}(i,:) = temp * this.nurbs_data_.control_points_(non_zero_id,1:3); end end %TODO computation of the Hassian matrix hassian = []; elseif any(strcmp(varargin, 'gradient')) for i = 1:size(xi,1) query_unit.query_protocol_ = {domain_patch, xi(i,:), 1}; this.basis_function_.query(query_unit); non_zero_id = query_unit.non_zero_id_; R = query_unit.evaluate_basis_{1}; position(i,:) = R * this.nurbs_data_.control_points_(non_zero_id,1:3); for dim_i = 1:geo_dim temp = query_unit.evaluate_basis_{2}(dim_i,:); gradient{dim_i}(i,:) = temp * this.nurbs_data_.control_points_(non_zero_id,1:3); end end elseif any(strcmp(varargin, 'position')) for i = 1:size(xi,1) query_unit.query_protocol_ = {domain_patch, xi(i,:), 0}; this.basis_function_.query(query_unit); non_zero_id = query_unit.non_zero_id_; R = query_unit.evaluate_basis_{1}; position(i,:) = R * this.nurbs_data_.control_points_(non_zero_id,1:3); end end end % Plot nurbs geometry function plotNurbs(this, varargin) geo_dim = this.nurbs_data_.getGeometryDimension(); if ~isempty(varargin) N = varargin{1}; else N = 11*ones(1, geo_dim); end switch geo_dim case 1 sample_pnt = linspace(0,1,N)'; position = this.evaluateNurbs(sample_pnt, 'position'); plot3(position(:,1), position(:,2), position(:,3), 'r-'); unique_knot = unique(this.nurbs_data_.knot_vectors{1}); position = this.evaluateNurbs(unique_knot, 'position'); plot3(position(:,1), position(:,2), position(:,3), 'k.'); case 2 unique_knot_vector{1} = unique(this.nurbs_data_.knot_vectors_{1}); unique_knot_vector{2} = unique(this.nurbs_data_.knot_vectors_{2}); this.plotSurfaceNurbs(N, unique_knot_vector); case 3 % Only surface nurbs are plotted for 3d case unique_knot_1 = unique(this.nurbs_data_.knot_vectors_{1}); unique_knot_2 = unique(this.nurbs_data_.knot_vectors_{2}); unique_knot_3 = unique(this.nurbs_data_.knot_vectors_{3}); % Plot face xi = 0 unique_knot_vector = {0, unique_knot_2, unique_knot_3}; this.plotSurfaceNurbs(N, unique_knot_vector); % Plot face xi = 1 unique_knot_vector = {1, unique_knot_2, unique_knot_3}; this.plotSurfaceNurbs(N, unique_knot_vector); % Plot face eta = 0 unique_knot_vector = {unique_knot_1, 0, unique_knot_3}; this.plotSurfaceNurbs(N, unique_knot_vector); % Plot face eta = 1 unique_knot_vector = {unique_knot_1, 1, unique_knot_3}; this.plotSurfaceNurbs(N, unique_knot_vector); % Plot face zeta = 0 unique_knot_vector = {unique_knot_1, unique_knot_2, 0}; this.plotSurfaceNurbs(N, unique_knot_vector); % Plot face zeta = 1 unique_knot_vector = {unique_knot_1, unique_knot_2, 1}; this.plotSurfaceNurbs(N, unique_knot_vect