CF round355 D题 构造 + 最小生成树

本文介绍了一种根据最小生成树(MST)边及非MST边的权值构造原始图的方法。通过将所有边按权值排序,并区分MST边与非MST边进行图构建,使用队列来辅助添加边,确保最终图的正确性和合理性。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

题意:

题目链接:http://codeforces.com/contest/606/problem/D
给出MST的边和非MST的边的权值,根据权值构造出原始的图。


思路:

先将所有的边按照权值排序,如果是MST的边就将点集内的点连接一个新的点,如果不是MST的边就在之前连接过的点集内添加一条边。
这里利用queue来保存可以用来构造的边,但是要注意设定一个上界2*m,否则会添加太多的边。


代码:

#include <bits/stdc++.h>
using namespace std;
const int MAXN = 1e5 + 10;

struct node {
    int v, t, id;
    bool operator < (const node &rhs) const {
        return v < rhs.v || (v == rhs.v && t > rhs.t);
    }
} e[MAXN];

int u[MAXN], v[MAXN];

int main() {
    //freopen("in.txt", "r", stdin);
    int n, m;
    scanf("%d%d", &n, &m);
    for (int i = 1; i <= m; i++) {
        scanf("%d%d", &e[i].v, &e[i].t);
        e[i].id = i;
    }
    sort(e + 1, e + 1 + m);
    int cnt = 1, num = 0;
    queue <long long> que;
    for (int i = 1; i <= m; i++) {
        int id = e[i].id;
        if (e[i].t == 1) {
            u[id] = cnt; v[id] = ++cnt;
            for (int j = cnt - 2; j >= 1 && num < 2 * m; j--, num++)
                que.push(j * 1000000LL + cnt);
        }
        else {
            if (que.empty()) {
                puts("-1");
                return 0;
            }
            long long tmp = que.front(); que.pop();
            u[id] = tmp / 1000000; v[id] = tmp % 1000000;
        }
    }
    for (int i = 1; i <= m; i++)
        printf("%d %d\n", u[i], v[i]);
    return 0;
}
# -*- coding: utf-8 -*- from abaqus import * from abaqusConstants import * from caeModules import * from odbAccess import * import math import csv import os import numpy as np from scipy.stats import qmc def run_fem_simulation(params): try: tool_position = params['tool_position'] tool_step = params['tool_step'] ae = params['ae'] R = params['R'] # ---------------几何参数定义---------------- ap = 10.0 # 轴向切深 ae = 3.0 # 径向切深 ft = 0.5 # 每齿进给量 min_mesh = 0.5 # 切削区网格尺寸 max_mesh = 3.0 # 非切削区网格尺寸 # ---------------材料参数定义---------------- length = 120.0 # 工件长度(X方向) width = 60.0 # 工件宽度(Z方向) depth = 5.0 # 工件厚度(Y方向) material_name = 'TC4' E = 113800.0 # 弹性模量 (MPa) nu = 0.33 # 泊松比 density = 4.43e-6 # 密度 (kg/mm³) # ---------------刀具参数定义---------------- R = 6.0 # 刀具半径 helix_angle1 = 38 # 螺旋角1,单位:度 helix_angle2 = 41 E_tool = 2.1e5 # MPa I = math.pi * R ** 4 / 4 # 圆截面惯性矩 # ---------------切削力系数------------------ kt1, kte1 = 2026.48, 7.16 # 切向力系数 kn1, kne1 = 972.23, 16.99 # 法向力系数 ka1, kae1 = 764.08, 0.51 # 轴向力系数 kt2, kte2 = 2294.84, 10.69 # 切向力系数 kn2, kne2 = 988.74, 18.20 # 法向力系数 ka2, kae2 = 1015.45, 1.39 # 轴向力系数 N = 4 # 刀具齿数 omega = 5000 # 主轴转速 (rpm) steps = 64 # 旋转步数 phis = math.pi / 2 - math.acos(1 - ae / R) phie = math.pi / 2 convergence_threshold = 0.1 # 收敛阈值 model_name = 'Model-%d-%d' % (tool_position, tool_step) if model_name in mdb.models.keys(): del mdb.models[model_name] myModel = mdb.Model(name=model_name) # ---------------几何建模--------------------- # Part创建 mySketch = myModel.ConstrainedSketch(name='sketch', sheetSize=200) mySketch.rectangle((0, 0), (length, depth)) myPart = myModel.Part(name='WorkPiece', dimensionality=THREE_D, type=DEFORMABLE_BODY) myPart.BaseSolidExtrude(sketch=mySketch, depth=width) # 分割工件(用于不同区域网格划分) datum = myPart.DatumPlaneByPrincipalPlane(principalPlane=XYPLANE, offset=width / 2) myPart.PartitionCellByDatumPlane( datumPlane=myPart.datums[datum.id], cells=myPart.cells.getSequenceFromMask(mask=('[#1 ]',), ) ) # ---------------网格划分--------------------- # 切削区域(上半部分)精细网格 region_edges = myPart.edges.getByBoundingBox( xMin=-1e-5, xMax=length + 1e-5, yMin=-1e-5, yMax=depth + 1e-5, zMin=width / 2 - 1e-5, zMax=width + 1e-5 ) myPart.seedEdgeBySize(edges=region_edges, size=min_mesh, constraint=FINER) # 非切削区域(下半部分)粗网格 other_edges = myPart.edges.getByBoundingBox( xMin=-1e-5, xMax=length + 1e-5, yMin=-1e-5, yMax=depth + 1e-5, zMin=-1e-5, zMax=width / 2 + 1e-5 ) myPart.seedEdgeBySize(edges=other_edges, size=max_mesh, constraint=FINER) myPart.generateMesh() # ---------------材料定义--------------------- myMaterial = myModel.Material(name=material_name) myMaterial.Density(table=((density,),)) myMaterial.Elastic(table=((E, nu),)) # 截面属性 mySection = myModel.HomogeneousSolidSection(name='Section', material=material_name) region = myPart.Set(cells=myPart.cells[:], name='AllCells') myPart.SectionAssignment(region=region, sectionName='Section') # ---------------装配体操作------------------- myAssembly = myModel.rootAssembly instance_name = 'WorkPiece-Instance' myInstance = myAssembly.Instance(name=instance_name, part=myPart, dependent=ON) # ---------------边界条件设置----------------- # 固定底面(Z=0) # 选择在 xz 平面上的所有节点 nodes = myInstance.nodes fixed_nodes = nodes.getByBoundingBox(xMin=-1e-05, xMax=length, yMin=-1e-05, yMax=depth, zMin=-1e-05, zMax=1e-05) fixed_node_set = myAssembly.Set(nodes=fixed_nodes, name='FixedNodeSet' ) myModel.DisplacementBC(name='BC-FIX', createStepName='Initial', region=fixed_node_set, u1=0.0, u2=0.0, u3=0.0 ) # ---------------分析步设置------------------- if tool_position > 1: myModel.StaticStep(name='Step-del', previous='Initial', maxNumInc=200, timePeriod=1e-05, initialInc=1e-10, minInc=1e-10, maxInc=1e-05 ) myModel.StaticStep(name='Step-1', previous='Step-del', nlgeom=ON, maxNumInc=200, timePeriod=0.1, initialInc=1e-05, minInc=1e-05, maxInc=1e-02 ) else: myModel.StaticStep(name='Step-1', previous='Initial', nlgeom=ON, maxNumInc=200, timePeriod=0.1, initialInc=1e-05, minInc=1e-05, maxInc=1e-02 ) ###---------------------------------------------创建被切削区域集合----------------------------------------------------- # 切削区节点坐标 xmin = -1e-05 xmax = length ymin = -1e-05 ymax = ae + 1e-05 zmin = width - ap + 1e-05 zmax = width + 1e-05 # 创建被切削节点集合 nodes = myInstance.nodes milling_nodes = nodes.getByBoundingBox(xMin=xmin, xMax=xmax, yMin=ymin, yMax=ymax, zMin=zmin, zMax=zmax) node_set_name = 'MillingNodes' myAssembly.Set(nodes=milling_nodes, name=node_set_name) print(f'Created node set [{node_set_name}] containing {len(milling_nodes)} nodes') # 创建被切削单元集合 elements = myInstance.elements milling_elements = elements.getByBoundingBox(xMin=xmin, xMax=xmax + min_mesh, yMin=ymin, yMax=ymax, zMin=zmin - min_mesh, zMax=zmax) element_set_name = 'MillingElements' myAssembly.Set(elements=milling_elements, name=element_set_name) print(f'Created element set [{element_set_name}] containing {len(milling_elements)} elements') #生死单元单元格集合 milling_cells = milling_elements.getByBoundingBox(xMin=- 1e-5, xMax=(tool_position - 1) * R + 1e-5, yMin=ymin, yMax=ymax, zMin=zmin - min_mesh, zMax=zmax) removed_name = f'remove-P{tool_position}' myAssembly.Set(elements=milling_cells, name=removed_name) print(f'Created removed feature set [{removed_name}] containing {len(milling_cells)} elements') #变形节点集合 deformation_nodes = nodes.getByBoundingBox(xMin=(tool_position - 1) * R - min_mesh/2, xMax=(tool_position - 1) * R + min_mesh/2, yMin=ae - min_mesh + min_mesh/2, yMax=ae + min_mesh - min_mesh/2, zMin=width - ap + tool_step * min_mesh - min_mesh/2, zMax=width - ap + tool_step * min_mesh + min_mesh/2) if deformation_nodes: deformation_set_name = f'DEFORMATION-P{tool_position}-M{tool_step}' myAssembly.Set(nodes=deformation_nodes, name=deformation_set_name) print(f'Created deformation set [{deformation_set_name}] with {len(deformation_nodes)} nodes') else: print(f'[Warning] No deformation nodes found for tool position {tool_position}') ###----------------------------------------------计算并施加载荷--------------------------------------------------------- convergence_threshold = 0.2 displacement = 0 previous_ae = 0 iteration_number = 1 if tool_position % 2 == 1: helix_angle = helix_angle1 else: helix_angle = helix_angle2 dphi = math.pi * 2 / steps if helix_angle != 0: db = R * dphi / math.tan(math.radians(helix_angle)) else: db = ap / steps steps_axial = int(ap / db + 1) milling_set = {} for node in milling_nodes: node_id = node.label node_x, node_y, node_z = node.coordinates milling_set[node_id] = (node_x, node_y, node_z) loaded_node_coords = [] # === 计算所有齿中底部微元(cnt4=0)的最小切出步 === dphi = math.pi * 2 / steps # 旋转步长 min_exit_step = steps # 初始化为最大值 for cnt3 in range(N): # 计算当前齿底部微元切出所需步数 phia_needed = math.pi / 2 - phis - cnt3 * 2 * math.pi / N # 调整phia_needed到0-2π范围内 while phia_needed < 0: phia_needed += 2 * math.pi while phia_needed >= 2 * math.pi: phia_needed -= 2 * math.pi # 计算切出步数 cnt1_exit = math.ceil(phia_needed / dphi) # 更新最小切出步 if cnt1_exit < min_exit_step: min_exit_step = cnt1_exit while abs(ae - previous_ae) > convergence_threshold: Fy_total = [0.0] * steps_axial zj_list = [width - ap + j * db for j in range(steps_axial)] # 每个cnt4对应的Z轴坐标 delta_y_current = [0.0] * steps_axial # 当前轮的刀具变形,初始化为0 loaded_node_coords = [] # 清空加载节点坐标 for cnt1 in range(steps): if cnt1 == (min_exit_step - 1) + tool_step: for cnt3 in range(N): if cnt3 in [0, 2]: # 齿编号0和2使用第一组系数 kt, kte = kt1, kte1 kn, kne = kn1, kne1 ka, kae = ka1, kae1 else: # 齿编号1和3使用第二组系数 kt, kte = kt2, kte2 kn, kne = kn2, kne2 ka, kae = ka2, kae2 for cnt4 in range(steps_axial): phia = phis + cnt1 * dphi + cnt3 * 2 * math.pi / N - cnt4 * db * math.tan( math.radians(helix_angle)) / R phia = phia % (2 * math.pi) if phis <= phia <= phie: h = ft * math.sin(phia) Ft = kt * h * db + kte * db Fn = kn * h * db + kne * db Fa = ka * h * db + kae * db else: Ft, Fn, Fa = 0, 0, 0 # 将切削力分解到直角坐标系 Fx = -Ft * math.cos(phia) - Fn * math.sin(phia) Fy = Ft * math.sin(phia) - Fn * math.cos(phia) Fz = Fa # 累加Fy Fy_total[cnt4] += Fy for i in range(int(length / R) - 1): # 计算刀具坐标系切削微元位置 micro_x = R * math.cos(phia) + (tool_position - 1) * R micro_y = R * math.sin(phia) - (depth - ae) - delta_y_current[cnt4] micro_z = width - ap + cnt4 * db # 找到最近的节点 min_distance = float('inf') nearest_node = None for node_id, (node_x, node_y, node_z) in milling_set.items(): distance = math.sqrt( (micro_x - node_x) ** 2 + (micro_y - node_y) ** 2 + (micro_z - node_z) ** 2) if distance < min_distance: min_distance = distance nearest_node = node_id # 创建载荷 if nearest_node is not None and (Fx != 0 or Fy != 0 or Fz != 0): node = milling_nodes.sequenceFromLabels((nearest_node,)) node_name = 'Load-' + str(cnt1) + '-' + str(cnt3) + '-' + str(cnt4) load_name = 'Force-' + str(cnt1) + '-' + str(cnt3) + '-' + str(cnt4) node_set = myAssembly.Set(nodes=node, name=node_name) myModel.ConcentratedForce(name=load_name, createStepName='Step-1', region=node_set, cf1=Fx, cf2=Fy, cf3=Fz) node_x, node_y, node_z = milling_set[nearest_node] loaded_node_coords.append((node_x, node_y, node_z)) # 创建参考点 adjusted_point = (micro_x, micro_y, micro_z) reference_point_name = 'RP-' + str(cnt1) + '-' + str(cnt3) + '-' + str(cnt4) reference_point = myAssembly.ReferencePoint(point=adjusted_point) # 刀具变形计算 delta_y_new = [] # 创建新列表存储本轮的变形计算结果 for j in range(steps_axial): zj = zj_list[j] delta = 0.0 for jp in range(steps_axial): zjp = zj_list[jp] if zj <= zjp: K = zjp ** 2 * (3 * zj - zjp) / (6 * E_tool * I) else: K = zj ** 2 * (3 * zjp - zj) / (6 * E_tool * I) delta += K * Fy_total[jp] delta_y_new.append(delta) # 更新delta_y_current用于下一次迭代(只在有计算值时更新) if delta_y_new: # 确保列表不为空 delta_y_current = delta_y_new ###------------------------------------------------------生死单元------------------------------------------------ if tool_position > 1: intname = 'Int-' + str(tool_position) stepname = 'Step-del' element_name = removed_name int_region = myAssembly.sets[element_name] myModel.ModelChange(name=intname, createStepName=stepname, region=int_region, regionType=ELEMENTS, activeInStep=False, includeStrain=False ) ###-------------------------------------------------------创建作业--------------------------------------------------- job_name = f'Model-{tool_position}-{tool_step}-{iteration_number}' myJob = mdb.Job(name=job_name, model=model_name, type=ANALYSIS, numCpus=24, numDomains=24) myJob.submit() mdb.jobs[job_name].waitForCompletion() ###-----------------------------------------------------打开ODB文件-------------------------------------------------- try: odb = openOdb(path=job_name + '.odb') step_name = 'Step-1' displacement_field = odb.steps[step_name].frames[-1].fieldOutputs['U'] instance_name = 'WORKPIECE-INSTANCE' instance = odb.rootAssembly.instances[instance_name] deformation_name = f'DEFORMATION-P{tool_position}-M{tool_step}' if deformation_name in odb.rootAssembly.nodeSets: region = odb.rootAssembly.nodeSets[deformation_name] displacement_subset = displacement_field.getSubset(region=region) displacement = displacement_subset.values[0].data[2] else: print(f"[Critical] Node set {deformation_name} missing in ODB!") displacement = 0 # 默认值或终止 print(f'ae={ae},displacement={displacement}') except Exception as e: print(f"Error reading ODB: {str(e)}") displacement = 0 finally: odb.close() pre_ae = ae ae -= abs(displacement) ae = (ae + previous_ae)/2 previous_ae = pre_ae if iteration_number > 10: break else: iteration_number += 1 final_displacement = displacement # 从ODB获取最终位移 final_ae = ae # 最终径向切深 print('final ae is:', final_ae) print('final displacement is:', final_displacement) return { 'displacement': final_displacement, 'final_ae': final_ae, 'status': 'success', 'message': 'Simulation completed successfully' } except Exception as e: return { 'status': 'error', 'message': str(e), 'displacement': None, 'final_ae': None } def generate_parameters(samples): """生成拉丁超立方采样参数组合""" param_ranges = { 'tool_position': (1, 19), 'tool_step': (1, 11) } # 拉丁超立方采样 sampler = qmc.LatinHypercube(d=len(param_ranges)) sample_data = sampler.random(n=samples) # 构造参数列表 params_list = [] for sample in qmc.scale(sample_data, [v[0] for v in param_ranges.values()], [v[1] for v in param_ranges.values()]): params = { 'tool_position': int(round(sample[0])), 'tool_step': int(round(sample[1])), 'ae': 3.0, # 固定值 'R': 5.0 # 固定值 } params_list.append(params) return params_list # 数据保存模块 class DataRecorder: def __init__(self, filename="simulation_data.csv"): self.filename = filename self._init_file() def _init_file(self): """初始化数据文件""" if not os.path.exists(self.filename): with open(self.filename, 'w', newline='') as f: writer = csv.writer(f) writer.writerow([ 'tool_position', 'tool_step', 'ae_input', 'R_input', 'displacement', 'final_ae', 'status' ]) def save_record(self, params, result): """保存单条记录""" record = { 'tool_position': params['tool_position'], 'tool_step': params['tool_step'], 'ae_input': params['ae'], 'R_input': params['R'], 'displacement': result.get('displacement', None), 'final_ae': result.get('final_ae', None), 'status': result['status'] } with open(self.filename, 'a', newline='') as f: writer = csv.DictWriter(f, fieldnames=record.keys()) writer.writerow(record) if __name__ == '__main__': # 初始化数据记录器 recorder = DataRecorder() # 生成参数组合 param_combinations = generate_parameters(samples=200) # 执行参数扫描 for idx, params in enumerate(param_combinations,1): print(f"\n 正在执行参数组合 {idx}/{len(param_combinations)}") print(f" 刀具位置: {params['tool_position']}") print(f" 刀具旋转步: {params['tool_step']}") print(f" 径向切深: {params['ae']}mm") print(f" 刀具半径: {params['R']}mm") # 执行仿真 result = run_fem_simulation(params) # 处理结果 if result['status'] == 'success': print(f"成功完成: 位移={result['displacement']}mm, 最终ae={result['final_ae']}mm") else: print(f"执行失败: {result['message']}") # 保存数据 recorder.save_record(params, result) print("\n所有参数组合执行完成!") print(f"数据文件已保存至: {os.path.abspath(recorder.filename)}") 以上是我的py代码,在abaqus中运行速度特别缓慢
最新发布
07-22
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值